/* 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
*
#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 "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()) {
-
- 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();
}
//________________________________________________________________________
//{
// // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
//
-// const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
+// AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance();
// TVector3 vecEmc ;
// TVector3 vecCpv ;
// if(cpv){
if(effPur>2 || effPur<0)
AliError(Form("Invalid Efficiency-Purity choice %d",effPur));
- if(ts->GetCpvIndex()<0)
- return 1 ; //no CPV cluster
+//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(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) ;
+// 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)))
// 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) ;
- //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 ;
const Int_t kSPECIES = AliPID::kSPECIESN ;
- AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
-
- Int_t nparticles = gime->RecParticles()->GetEntriesFast() ;
+ Int_t nparticles = fRecParticles->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)
- 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. ;
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() ;
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).
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);
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<AliPHOSTrackSegment *>(fTrackSegments->At(rp->GetPHOSTSIndex()));
+ AliPHOSEmcRecPoint * erp = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(ts->GetEmcIndex()));
TVector3 pos ;
- geom->GetGlobal(erp, pos) ;
+ fGeom->GetGlobalPHOS(erp, pos) ;
rp->SetPos(pos);
index++ ;
}
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() ;
}
{
// 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 += " " ;
}
//____________________________________________________________________________
-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
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff