-#include "AliReaderESD.h"
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
//____________________________________________________________________
//////////////////////////////////////////////////////////////////////
// //
// //
//////////////////////////////////////////////////////////////////////
-#include <TPDGCode.h>
-#include <TString.h>
-#include <TObjString.h>
-#include <TTree.h>
#include <TFile.h>
-#include <TKey.h>
-#include <TParticle.h>
+#include <TGliteXmlEventlist.h>
#include <TH1.h>
+#include <TPDGCode.h>
+#include <TParticle.h>
+#include <TString.h>
-#include <AliRun.h>
-#include <AliRunLoader.h>
-#include <AliStack.h>
-#include <AliESDtrack.h>
-#include <AliKalmanTrack.h>
-#include <AliESD.h>
-
-#include "AliAnalysis.h"
-#include "AliAODRun.h"
#include "AliAOD.h"
#include "AliAODParticle.h"
-#include "AliAODParticleCut.h"
-#include "AliTrackPoints.h"
#include "AliClusterMap.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliLog.h"
+#include "AliReaderESD.h"
+#include "AliRunLoader.h"
+#include "AliStack.h"
ClassImp(AliReaderESD)
fdR(0.0),
fClusterMap(kFALSE),
fITSTrackPoints(kFALSE),
+ fITSTrackPointsType(AliTrackPoints::kITS),
fMustTPC(kFALSE),
+ fReadCentralBarrel(kTRUE),
+ fReadMuon(kFALSE),
+ fReadPHOS(kFALSE),
fNTPCClustMin(0),
- fNTPCClustMax(150),
+ fNTPCClustMax(1500),
fTPCChi2PerClustMin(0.0),
fTPCChi2PerClustMax(10e5),
fChi2Min(0.0),
{
//cosntructor
- if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
+ if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
}
/********************************************************************/
fdR(0.0),
fClusterMap(kFALSE),
fITSTrackPoints(kFALSE),
+ fITSTrackPointsType(AliTrackPoints::kITS),
fMustTPC(kFALSE),
+ fReadCentralBarrel(kTRUE),
+ fReadMuon(kFALSE),
+ fReadPHOS(kFALSE),
fNTPCClustMin(0),
fNTPCClustMax(150),
fTPCChi2PerClustMin(0.0),
fTPCC44Max(10e5)
{
//cosntructor
- if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
+ if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
}
/********************************************************************/
AliReaderESD::~AliReaderESD()
{
//desctructor
- delete fRunLoader;
- delete fKeyIterator;
- delete fFile;
+ delete fRunLoader;
+ delete fKeyIterator;
+ delete fFile;
}
+
/**********************************************************/
Int_t AliReaderESD::ReadNext()
{
//reads next event from fFile
//fRunLoader is for reading Kine
- if (AliVAODParticle::GetDebug())
- Info("ReadNext","Entered");
+ AliDebug(1,"Entered");
if (fEventSim == 0x0) fEventSim = new AliAOD();
if (fEventRec == 0x0) fEventRec = new AliAOD();
fEventRec->Reset();
do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
- {
- if (fFile == 0x0)
- {
- fFile = OpenFile(fCurrentDir);//rl is opened here
- if (fFile == 0x0)
- {
- Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
- fCurrentDir++;
- continue;
- }
- fCurrentEvent = 0;
- fKeyIterator = new TIter(fFile->GetListOfKeys());
-// fFile->Dump();
-// fFile->GetListOfKeys()->Print();
- }
- TKey* key = (TKey*)fKeyIterator->Next();
- if (key == 0x0)
- {
- if (AliVAODParticle::GetDebug() > 2 )
- {
- Info("ReadNext","No more keys.");
- }
- fCurrentDir++;
- delete fKeyIterator;
- fKeyIterator = 0x0;
- delete fFile;//we have to assume there is no more ESD objects in the fFile
- fFile = 0x0;
- delete fRunLoader;
- fRunLoader = 0x0;
- continue;
- }
- //try to read
-
-
-// TObject* esdobj = key->ReadObj();
-// if (esdobj == 0x0)
-// {
-// if (AliVAODParticle::GetDebug() > 2 )
-// {
-// Info("ReadNext","Key read NULL. Key Name is %s",key->GetName());
-// key->Dump();
-// }
-// continue;
-// }
-// esdobj->Dump();
-// AliESD* esd = dynamic_cast<AliESD*>(esdobj);
-
+ {
+ if (fFile == 0x0)
+ {
+ fFile = OpenFile(fCurrentDir);//rl is opened here
+ if (fFile == 0x0)
+ {
+ Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
+ fCurrentDir++;
+ continue;
+ }
+ fCurrentEvent = 0;
+ }
TString esdname = "ESD";
esdname+=fCurrentEvent;
AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
if (esd == 0x0)
{
-// if (AliVAODParticle::GetDebug() > 2 )
-// {
-// Info("ReadNext","This key is not an AliESD object %s",key->GetName());
-// }
- if (AliVAODParticle::GetDebug() > 2 )
- {
- Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
- }
+ AliDebug(3,Form("Can not find AliESD object named %s",esdname.Data()));
fCurrentDir++;
- delete fKeyIterator;
- fKeyIterator = 0x0;
delete fFile;//we have to assume there is no more ESD objects in the fFile
fFile = 0x0;
delete fRunLoader;
fRunLoader = 0x0;
continue;
}
-
- ReadESD(esd);
+ ReadESD(esd);
- fCurrentEvent++;
- fNEventsRead++;
- delete esd;
- return 0;//success -> read one event
- }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
+ fCurrentEvent++;
+ fNEventsRead++;
+ delete esd;
+ return 0;//success -> read one event
+ }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
return 1; //no more directories to read
}
-/**********************************************************/
+/**********************************************************/
Int_t AliReaderESD::ReadESD(AliESD* esd)
+{
+//Reads esd data
+ if (esd == 0x0)
+ {
+ Error("ReadESD","ESD is NULL");
+ return 1;
+ }
+
+ // seperate each method
+ if (fReadCentralBarrel) ReadESDCentral(esd);
+
+ if (fReadMuon) ReadESDMuon(esd);
+
+ if (fReadPHOS) ReadESDPHOS(esd);
+
+ return 1;
+}
+
+/**********************************************************/
+Int_t AliReaderESD::ReadESDCentral(AliESD* esd)
{
//****** Tentative particle type "concentrations"
- static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
+ static const Double_t kConcentr[5]={0.05, 0., 0.85, 0.10, 0.05};
Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
Double_t w[kNSpecies];
Double_t pos[3];//position
Double_t vertexpos[3];//vertex position
//Reads one ESD
- if (esd == 0x0)
- {
- Error("ReadESD","ESD is NULL");
- return 1;
- }
TDatabasePDG* pdgdb = TDatabasePDG::Instance();
if (pdgdb == 0x0)
return 1;
}
- Float_t mf = esd->GetMagneticField();
+ Float_t mf = esd->GetMagneticField()/10.; //AliESD::GetMagnField returns mf in kG
if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) )
{
Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
- return 1;
+
}
if (fITSTrackPoints)
{
- Info("ReadESD","Magnetic Field is %f",mf/10.);
- AliKalmanTrack::SetMagneticField(mf/10.);
+ Info("ReadESD","Magnetic Field is %f",mf);
+ //AliKalmanTrack::SetMagneticField(mf);
}
AliStack* stack = 0x0;
vertex->GetXYZ(vertexpos);
}
- if (AliVAODParticle::GetDebug() > 0)
- {
- Info("ReadESD","Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]);
- }
+ AliDebug(1,Form("Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]));
Info("ReadESD","Reading Event %d",fCurrentEvent);
//if (esdtrack->HasVertexParameters() == kFALSE)
if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
{
- if (AliVAODParticle::GetDebug() > 2)
- Info("ReadNext","Particle skipped: Data at vertex not available.");
+ AliDebug(3,Form("Particle skipped: Data at vertex not available."));
continue;
}
{
if ((esdtrack->GetStatus() & AliESDtrack::kTPCin) == kFALSE)
{
- if (AliVAODParticle::GetDebug() > 2)
- Info("ReadNext","Particle skipped: Was not reconstructed in TPC.");
+ AliDebug(3,"Particle skipped: Was not reconstructed in TPC.");
continue;
}
}
if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE)
{
- if (AliVAODParticle::GetDebug() > 2)
- Info("ReadNext","Particle skipped: PID BIT is not set.");
+ AliDebug(3,"Particle skipped: PID BIT is not set.");
continue;
}
- Double_t extx;
+ Double_t alpha,extx;
Double_t extp[5];
- esdtrack->GetConstrainedExternalParameters(extx,extp);
+ esdtrack->GetConstrainedExternalParameters(alpha,extx,extp);
if (extp[4] == 0.0)
{
- if (AliVAODParticle::GetDebug() > 2)
- Info("ReadNext","Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
+ AliDebug(3,"Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
continue;
}
esdtrack->GetESDpid(pidtable);
{
if(Rejected(p->GetPdgCode()))
{
- if ( AliVAODParticle::GetDebug() > 5 )
- Info("ReadNext","Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode());
+ AliDebug(6,Form("Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode()));
continue; //check if we are intersted with particles of this type
}
}
//Here we apply Bayes' formula
Double_t rc=0.;
- for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
+ for (Int_t s=0; s<AliPID::kSPECIES; s++) rc+=kConcentr[s]*pidtable[s];
if (rc==0.0)
{
- if (AliVAODParticle::GetDebug() > 2)
- Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
+ AliDebug(3,"Particle rejected since total bayessian PID probab. is zero.");
continue;
}
- for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
+ for (Int_t s=0; s<AliPID::kSPECIES; s++) w[s]=kConcentr[s]*pidtable[s]/rc;
- if (AliVAODParticle::GetDebug() > 4)
+ if (AliDebugLevel() > 4)
{
Info("ReadNext","###########################################################################");
Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
msg+=")";
}
Info("ReadNext","%s",msg.Data());
- }//if (AliVAODParticle::GetDebug()>4)
+ }//if (AliDebugLevel()>4)
AliTrackPoints* tpts = 0x0;
if (fNTrackPoints > 0)
{
- tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf,fdR);
+ tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf*10.0,fdR);
// tpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
}
AliTrackPoints* itstpts = 0x0;
if (fITSTrackPoints)
{
- itstpts = new AliTrackPoints(AliTrackPoints::kITS,esdtrack);
+ itstpts = new AliTrackPoints(fITSTrackPointsType,esdtrack,mf*10.0);
// itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
}
//find the most probable PID
Int_t spec = 0;
Float_t maxprob = w[0];
- for (Int_t s=1; s<AliESDtrack::kSPECIES; s++)
+ for (Int_t s=1; s<AliPID::kSPECIES; s++)
{
if (w[s]>maxprob)
{
Float_t pp = w[s];
if (pp == 0.0)
{
- if ( AliVAODParticle::GetDebug() > 5 )
- Info("ReadNext","Probability of being PID %d is zero. Continuing.",pdgcode);
+ AliDebug(6,Form("Probability of being PID %d is zero. Continuing.",pdgcode));
continue;
}
if(Rejected(pdgcode))
{
- if ( AliVAODParticle::GetDebug() > 5 )
- Info("ReadNext","PID (%d) did not pass the cut.",pdgcode);
+ AliDebug(6,Form("PID (%d) did not pass the cut.",pdgcode));
continue; //check if we are intersted with particles of this type
}
if(Rejected(track))//check if meets all criteria of any of our cuts
//if it does not delete it and take next good track
{
- if ( AliVAODParticle::GetDebug() > 4 )
- Info("ReadNext","Track did not pass the cut");
+ AliDebug(5,"Track did not pass the cut");
delete track;
continue;
}
if (particle) fEventSim->AddParticle(particle);
keeppart = kTRUE;
- if (AliVAODParticle::GetDebug() > 4 )
+ if (AliDebugLevel() > 4 )
{
Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
track->Print();
}
else
{
- if (particle->P() < 0.00001)
+ if ( fReadSim && stack )
{
- Info("ReadNext","###################################");
- Info("ReadNext","###################################");
- Info("ReadNext","Track Label %d",esdtrack->GetLabel());
- TParticle *p = stack->Particle(esdtrack->GetLabel());
- Info("ReadNext","");
- p->Print();
- Info("ReadNext","");
- particle->Print();
- }
+ if (particle->P() < 0.00001)
+ {
+ Info("ReadNext","###################################");
+ Info("ReadNext","###################################");
+ Info("ReadNext","Track Label %d",esdtrack->GetLabel());
+ TParticle *p = stack->Particle(esdtrack->GetLabel());
+ Info("ReadNext","");
+ p->Print();
+ Info("ReadNext","");
+ particle->Print();
+ }
+ }
}
}//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
return 0;
}
+/**********************************************************/
+Int_t AliReaderESD::ReadESDMuon(AliESD* esd)
+{
+ // Reads the muon tracks from the ESD
+ Double_t vertexpos[3];//vertex position, assuming no secondary decay
+
+ const AliESDVertex* vertex = esd->GetVertex();
+
+ if (vertex == 0x0) {
+ Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
+ vertexpos[0] = 0.0;
+ vertexpos[1] = 0.0;
+ vertexpos[2] = 0.0;
+ } else {
+ vertex->GetXYZ(vertexpos);
+ }
+
+ Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
+
+ AliDebug(1,Form("Reading Event %d \nFound %d tracks.",fCurrentEvent,nTracks));
+
+ // settings
+ Float_t chi2Cut = 100.;
+ Float_t ptCutMin = 1.;
+ Float_t ptCutMax = 10000.;
+ Float_t muonMass = 0.105658389;
+ Int_t pdgcode = -13;
+ Double_t thetaX, thetaY, pYZ;
+ Double_t pxRec1, pyRec1, pzRec1, e1;
+ Int_t charge;
+
+ Int_t ntrackhits;
+ Double_t fitfmin;
+
+ TLorentzVector fV1;
+ fEventRec->Reset();
+ for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
+
+ AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
+
+ thetaX = muonTrack->GetThetaX();
+ thetaY = muonTrack->GetThetaY();
+
+ pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
+ pzRec1 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX));
+ pxRec1 = pzRec1 * TMath::Tan(thetaX);
+ pyRec1 = pzRec1 * TMath::Tan(thetaY);
+ charge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
+ e1 = TMath::Sqrt(muonMass * muonMass + pxRec1 * pxRec1 + pyRec1 * pyRec1 + pzRec1 * pzRec1);
+ fV1.SetPxPyPzE(pxRec1, pyRec1, pzRec1, e1);
+
+ ntrackhits = muonTrack->GetNHit();
+ fitfmin = muonTrack->GetChi2();
+
+ // transverse momentum
+ Float_t pt1 = fV1.Pt();
+
+ // chi2 per d.o.f.
+ Float_t ch1 = fitfmin / (2.0 * ntrackhits - 5);
+
+ if ((ch1 < chi2Cut) && (pt1 > ptCutMin) && (pt1 < ptCutMax)) {
+ AliAODParticle* track = new AliAODParticle(pdgcode*charge,1,iTrack,
+ pxRec1, pyRec1,pzRec1, e1,
+ vertexpos[0], vertexpos[1], vertexpos[2], 0.);
+ fEventRec->AddParticle(track);
+ }
+
+ }
+ fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
+ return 0;
+}
+
/**********************************************************/
void AliReaderESD::Rewind()
{
//rewinds reading
- delete fKeyIterator;
+ // delete fKeyIterator;
delete fFile;
fFile = 0x0;
- fKeyIterator = 0x0;
+ // fKeyIterator = 0x0;
delete fRunLoader;
fRunLoader = 0x0;
fCurrentDir = 0;
fNEventsRead = 0;
+ if (fEventList) fEventList->Reset();
if (fTrackCounter) fTrackCounter->Reset();
}
/**********************************************************/
TFile* AliReaderESD::OpenFile(Int_t n)
{
-//opens fFile with kine tree
-
+//opens fFile with tree
+ if (fEventList)
+ {
+ if (fCurrentDir > n)
+ {
+ fEventList->Reset();
+ fCurrentDir = 0;
+ }
+
+ while (fCurrentDir < n)
+ {
+ fEventList->Next();
+ fCurrentDir++;
+ }
+ fEventList->Next();
+ }
+
const TString& dirname = GetDirName(n);
if (dirname == "")
{
Error("OpenFiles","Can not get directory name");
return 0x0;
}
- TString filename = dirname +"/"+ fESDFileName;
+ TString filename;
+
+ if (fEventList)
+ {
+ filename = fEventList->GetURL(fESDFileName);
+ }
+ else
+ {
+ filename = dirname +"/"+ fESDFileName;
+ }
+ Info("OpenFile","%s ==> %s",fESDFileName.Data(),filename.Data());
TFile *ret = TFile::Open(filename.Data());
if ( ret == 0x0)
if (fReadSim )
{
- fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
+ TString gafilename;
+ if (fEventList)
+ {
+ gafilename = fEventList->GetURL(fGAlFileName);
+ }
+ else
+ {
+ gafilename = dirname +"/"+ fGAlFileName;
+ }
+ Info("OpenFile","%s ==> %s",fGAlFileName.Data(),gafilename.Data());
+
+ fRunLoader = AliRunLoader::Open(gafilename);
+
if (fRunLoader == 0x0)
{
Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
}
fRunLoader->LoadHeader();
+
+ if (fEventList)
+ {
+ TString kinefilename = fEventList->GetURL("Kinematics.root");
+ fRunLoader->SetKineFileName(kinefilename);
+ Info("OpenFile","%s ==> %s","Kinematics.root",kinefilename.Data());
+ }
+
if (fRunLoader->LoadKinematics())
{
Error("Next","Error occured while loading kinematics.");
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff