X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ANALYSIS%2FAliReaderESD.cxx;h=b6aa2dfe1fb79eb9ae14da4afe15a509ec422d84;hb=abda854188ebf50082975185957c100147a929a1;hp=cc02a73f9082ed30ad84b34bb3a14ca77a15e450;hpb=afa8b37b24ee004cb6b499d65baaea4476318cd1;p=u%2Fmrichter%2FAliRoot.git diff --git a/ANALYSIS/AliReaderESD.cxx b/ANALYSIS/AliReaderESD.cxx index cc02a73f908..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,28 +26,22 @@ // // ////////////////////////////////////////////////////////////////////// -#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) @@ -47,8 +57,14 @@ AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename) fNTrackPoints(0), 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), @@ -76,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."); } /********************************************************************/ @@ -94,6 +110,12 @@ AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char fNTrackPoints(0), 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), @@ -122,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."); } /********************************************************************/ @@ -130,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(); @@ -150,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]; @@ -243,11 +240,6 @@ 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) @@ -256,14 +248,20 @@ 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) ) + 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); + //AliKalmanTrack::SetMagneticField(mf); + } + AliStack* stack = 0x0; if (fReadSim && fRunLoader) { @@ -284,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); @@ -305,26 +300,32 @@ 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; } + if (fMustTPC) + { + if ((esdtrack->GetStatus() & AliESDtrack::kTPCin) == kFALSE) + { + 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); @@ -348,12 +349,15 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) Error("ReadNext","Can not find track with such label."); continue; } - if(Pass(p->GetPdgCode())) + + if (fCheckParticlePID) { - if ( AliVAODParticle::GetDebug() > 5 ) - Info("ReadNext","Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode()); - continue; //check if we are intersted with particles of this type - } + if(Rejected(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 + } + } // if(p->GetPdgCode()<0) charge = -1; particle = new AliAODParticle(*p,i); @@ -363,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]); @@ -392,12 +395,20 @@ 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(fITSTrackPointsType,esdtrack,mf*10.0); +// itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]); } AliClusterMap* cmap = 0x0; @@ -419,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) { @@ -437,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(Pass(pdgcode)) + 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 } @@ -463,11 +472,10 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]); } - if(Pass(track))//check if meets all criteria of any of our cuts + 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; } @@ -475,7 +483,12 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) //Single Particle cuts on cluster map and track points rather do not have sense if (tpts) { - track->SetTrackPoints(tpts); + track->SetTPCTrackPoints(tpts); + } + + if (itstpts) + { + track->SetITSTrackPoints(itstpts); } if (cmap) @@ -487,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(); @@ -500,8 +513,26 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) { delete particle;//particle was not stored in event delete tpts; + delete itstpts; delete cmap; } + else + { + if ( fReadSim && stack ) + { + 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;iGetNumberOfParticles(), fEventSim->GetNumberOfParticles(), fNEventsRead,fCurrentEvent,fCurrentDir); fTrackCounter->Fill(fEventRec->GetNumberOfParticles()); + + /******************************************************/ + /****** Setting glevet properties *************/ + /******************************************************/ + + if (fEventRec->GetNumberOfParticles() > 0) + { + fEventRec->SetPrimaryVertex(vertexpos[0],vertexpos[1],vertexpos[2]); + } + 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; } @@ -517,29 +629,55 @@ Int_t AliReaderESD::ReadESD(AliESD* esd) 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) @@ -555,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()); @@ -564,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.");