X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ANALYSIS%2FAliReaderESD.cxx;h=b6aa2dfe1fb79eb9ae14da4afe15a509ec422d84;hb=ec6465fe79f359cc2601100cb4a851a547fa08a8;hp=17e09a737d065511382c807bfc2a584a3a7cd28f;hpb=f4b459f7a53843ddce7b2844b4fb2304e3ab9367;p=u%2Fmrichter%2FAliRoot.git diff --git a/ANALYSIS/AliReaderESD.cxx b/ANALYSIS/AliReaderESD.cxx index 17e09a737d0..b6aa2dfe1fb 100644 --- a/ANALYSIS/AliReaderESD.cxx +++ b/ANALYSIS/AliReaderESD.cxx @@ -1,4 +1,20 @@ -#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$ */ + //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// // // @@ -10,29 +26,22 @@ // // ////////////////////////////////////////////////////////////////////// -#include -#include -#include -#include #include -#include -#include +#include #include +#include +#include +#include -#include -#include -#include -#include -#include -#include - -#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) @@ -49,9 +58,13 @@ AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename) 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), @@ -79,7 +92,7 @@ AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename) { //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."); } /********************************************************************/ @@ -98,7 +111,11 @@ AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char 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), @@ -127,7 +144,7 @@ AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char 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."); } /********************************************************************/ @@ -135,18 +152,18 @@ AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char 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(); @@ -155,92 +172,67 @@ Int_t AliReaderESD::ReadNext() 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(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(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]; @@ -248,12 +240,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) 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) { @@ -261,18 +248,18 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) 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; @@ -295,10 +282,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) 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); @@ -316,8 +300,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) //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; } @@ -325,27 +308,24 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) { 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); @@ -374,8 +354,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) { 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 } } @@ -388,17 +367,16 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) //Here we apply Bayes' formula Double_t rc=0.; - for (Int_t s=0; s 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 4) + if (AliDebugLevel() > 4) { Info("ReadNext","###########################################################################"); Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]); @@ -417,19 +395,19 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) 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]); } @@ -452,7 +430,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) //find the most probable PID Int_t spec = 0; Float_t maxprob = w[0]; - for (Int_t s=1; smaxprob) { @@ -470,15 +448,13 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) 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 } @@ -499,8 +475,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) 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; } @@ -525,7 +500,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) 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(); @@ -577,34 +552,132 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) 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) @@ -620,7 +693,19 @@ TFile* AliReaderESD::OpenFile(Int_t n) 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()); @@ -629,6 +714,14 @@ TFile* AliReaderESD::OpenFile(Int_t n) } 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.");