/* 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
*
#include "TPrincipal.h"
#include "TFile.h"
#include "TSystem.h"
-#include "TVector3.h"
// --- AliRoot header files ---
//#include "AliLog.h"
#include "AliPHOS.h"
#include "AliPHOSPIDv1.h"
-#include "AliPHOSGetter.h"
-#include "AliESD.h"
+#include "AliESDEvent.h"
#include "AliESDVertex.h"
+#include "AliPHOSTrackSegment.h"
+#include "AliPHOSEmcRecPoint.h"
+#include "AliPHOSRecParticle.h"
ClassImp( AliPHOSPIDv1)
//____________________________________________________________________________
AliPHOSPIDv1::AliPHOSPIDv1() :
+ AliPHOSPID(),
fBayesian(kFALSE),
fDefaultInit(kFALSE),
fWrite(kFALSE),
- fNEvent(0),
fFileNamePrincipalPhoton(),
fFileNamePrincipalPi0(),
fFileNameParameters(),
fX(0),
fPPhoton(0),
fPPi0(0),
- fRecParticlesInRun(0),
fParameters(0),
+ fVtx(0.),
fTFphoton(0),
fTFpiong(0),
fTFkaong(0),
fBayesian(kFALSE),
fDefaultInit(kFALSE),
fWrite(kFALSE),
- fNEvent(0),
fFileNamePrincipalPhoton(),
fFileNamePrincipalPi0(),
fFileNameParameters(),
fX(0),
fPPhoton(0),
fPPi0(0),
- fRecParticlesInRun(0),
fParameters(0),
+ fVtx(0.),
fTFphoton(0),
fTFpiong(0),
fTFkaong(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),
- fNEvent(0),
fFileNamePrincipalPhoton(),
fFileNamePrincipalPi0(),
fFileNameParameters(),
fX(0),
fPPhoton(0),
fPPi0(0),
- fRecParticlesInRun(0),
fParameters(0),
+ fVtx(0.),
fTFphoton(0),
fTFpiong(0),
fTFkaong(0),
//ctor with the indication on where to look for the track segments
InitParameters() ;
- Init() ;
fDefaultInit = kFALSE ;
}
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
}
//________________________________________________________________________
-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");
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()) {
-
- GetVertex() ;
- 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();
}
//________________________________________________________________________
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
{
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 ;
-//}
//____________________________________________________________________________
Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSTrackSegment * ts, Int_t effPur, Float_t e) 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
// 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) ;
- TMatrixF dummy ;
-
- emc->GetGlobalPosition(dir, dummy) ;
-
TVector3 local ;
emc->GetLocalPosition(local) ;
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() ;
+ 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{
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() ;
+ depthx = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
+ depthz = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
}
- 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 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) ;
- //Correct for the vertex position
- dir = dir - fVtx ;
+ dir-=fVtx ;
dir.SetMag(1.) ;
return dir ;
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 ;
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)
//// track = ts->GetTrackIndex() ; //TPC tracks ?
if (!emc) {
- AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ;
+ AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ;
}
for(index = 0 ; index < nparticles ; index ++) {
- AliPHOSRecParticle * recpar = gime->RecParticle(index) ;
+ AliPHOSRecParticle * recpar = static_cast<AliPHOSRecParticle *>(fRecParticles->At(index));
//Conversion electron?
{
// 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 "<<index<<endl;
- new( (*recParticles)[index] ) AliPHOSRecParticle() ;
- rp = (AliPHOSRecParticle *)recParticles->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() ;
// 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)){
rp->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(fVtx.X(),fVtx.Y(),fVtx.Z(),0);
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<AliPHOSTrackSegment *>(fTrackSegments->At(rp->GetPHOSTSIndex()));
+ AliPHOSEmcRecPoint * erp = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(ts1->GetEmcIndex()));
TVector3 pos ;
- geom->GetGlobal(erp, pos) ;
+ fGeom->GetGlobalPHOS(erp, pos) ;
rp->SetPos(pos);
index++ ;
}
{
// 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 += " " ;
&(*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);
}
(*fParameters)(p,i) = 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){
+ 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.) ;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff