X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=JETAN%2FAliJetAODReader.cxx;h=2bdf7590f4dc75cf1147c15cd43caaae33d03d53;hp=597cf7defa23629ef605eedb68a34b263f67b540;hb=4bd66a92319fe2800eedcf7180ae924178dc2071;hpb=9f16612007431eadab205efc5883154499d2d81d diff --git a/JETAN/AliJetAODReader.cxx b/JETAN/AliJetAODReader.cxx index 597cf7defa2..2bdf7590f4d 100644 --- a/JETAN/AliJetAODReader.cxx +++ b/JETAN/AliJetAODReader.cxx @@ -20,6 +20,11 @@ // is executed after the ESD filter task, in order to read its output // // Author: Davide Perrino +// +// **February 2011 +// implemented standard geometry (AliEMCALGeometry) instead of dummy one (AliJetDummyGeo) +// moved geometry definition in AliJetReader +// marco.bregant@subatech.in2p3.fr //------------------------------------------------------------------------- @@ -31,18 +36,22 @@ #include #include #include +#include #include "AliJetAODReader.h" #include "AliJetAODReaderHeader.h" #include "AliAODEvent.h" #include "AliAODTrack.h" #include "AliAODMCParticle.h" -#include "AliJetDummyGeo.h" +#include "AliEMCALGeometry.h" #include "AliJetAODFillUnitArrayTracks.h" #include "AliJetAODFillUnitArrayEMCalDigits.h" #include "AliJetHadronCorrectionv1.h" #include "AliJetUnitArray.h" +#include "AliOADBContainer.h" +using std::cout; +using std::endl; ClassImp(AliJetAODReader) AliJetAODReader::AliJetAODReader(): @@ -51,7 +60,6 @@ AliJetAODReader::AliJetAODReader(): fRef(new TRefArray), fDebug(0), fOpt(0), - fGeom(0), fHadCorr(0x0), fTpcGrid(0x0), fEmcalGrid(0x0), @@ -119,9 +127,7 @@ void AliJetAODReader::OpenInputFiles() if (a>=naod) continue; if (strstr(name,pattern)){ - char path[256]; - sprintf(path,"%s/%s/aod.root",dirName,name); - fChain->AddFile(path); + fChain->AddFile(Form("%s/%s/aod.root",dirName,name)); a++; } } @@ -267,40 +273,118 @@ Bool_t AliJetAODReader::FillMomentumArray() return kFALSE; } - if(((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC()){ - return FillMomentumArrayMC(); - } - - // get number of tracks in event (for the loop) - Int_t nt = fAOD->GetNTracks(); - if(fDebug>0)printf("AOD tracks: %5d \t", nt); - - // temporary storage of signal and pt cut flag - Int_t* sflag = new Int_t[nt]; - Int_t* cflag = new Int_t[nt]; - // get cuts set by user Float_t ptMin = fReaderHeader->GetPtCut(); Float_t etaMin = fReaderHeader->GetFiducialEtaMin(); Float_t etaMax = fReaderHeader->GetFiducialEtaMax(); UInt_t filterMask = ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask(); - //loop over tracks + // ----- number of tracks ----- + Int_t nTracksStd = 0; + Short_t mcReaderFlag = ((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC(); + TClonesArray *mcArray = 0x0; + // check if we have to read from MC branch + if (((AliJetAODReaderHeader*)fReaderHeader)->GetReadStdBranch()) { + if(mcReaderFlag) { + mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName()); + if(!mcArray){ + Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__); + } + nTracksStd = mcArray->GetEntriesFast(); + } + else { + nTracksStd = fAOD->GetNTracks(); + // printf("no. of standard tracks: %i\n", nTracksStd); + } + } + Int_t nTracksNonStd = 0; + TClonesArray *nonStdTracks = 0x0; + if (((AliJetAODReaderHeader*)fReaderHeader)->GetReadNonStdBranch()) { + nonStdTracks = + (TClonesArray*) fAOD->FindListObject(((AliJetAODReaderHeader*)fReaderHeader)->GetNonStdBranch()); + if (nonStdTracks) + nTracksNonStd = nonStdTracks->GetEntries(); + // printf("no. of non-standard tracks: %i\n", nTracksNonStd); + } + Int_t nTracks = nTracksStd + nTracksNonStd; + + // temporary storage of signal and pt cut flag + Int_t* sflag = new Int_t[nTracks]; + Int_t* cflag = new Int_t[nTracks]; + + // loop over tracks Int_t aodTrack = 0; Float_t pt, eta; TVector3 p3; - for (Int_t it = 0; it < nt; it++) { - AliAODTrack *track = fAOD->GetTrack(it); - UInt_t status = track->GetStatus(); - + // ----- looping over standard branch ----- + if (mcArray) { + for (Int_t iTrack = 0; iTrack < nTracksStd; iTrack++) { + AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(iTrack); + if(!track->IsPhysicalPrimary())continue; + + p3.SetXYZ(track->Px(),track->Py(),track->Pz()); + eta = p3.Eta(); + if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut + if(!AcceptAODMCParticle(track,mcReaderFlag))continue; + pt = p3.Pt(); + + if (aodTrack == 0){ + fRef->Delete(); // make sure to delete before placement new... + new(fRef) TRefArray(TProcessID::GetProcessWithUID(track)); + } + new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag()); + sflag[aodTrack] = 1; + cflag[aodTrack] = ( pt > ptMin ) ? 1: 0; + fRef->Add(track); + aodTrack++; + } + } + else { + for (Int_t iTrack = 0; iTrack < nTracksStd; iTrack++) { + AliAODTrack *track = fAOD->GetTrack(iTrack); + UInt_t status = track->GetStatus(); + + Double_t mom[3] = {track->Px(),track->Py(),track->Pz()}; + p3.SetXYZ(mom[0],mom[1],mom[2]); + pt = p3.Pt(); + eta = p3.Eta(); + if (status == 0) continue; + if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue; + if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut + + if (aodTrack == 0){ + fRef->Delete(); // make sure to delete before placement new... + new(fRef) TRefArray(TProcessID::GetProcessWithUID(track)); + } + new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag()); + sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0; + cflag[aodTrack] = ( pt > ptMin ) ? 1: 0; + aodTrack++; + fRef->Add(track); + } + } + + // ----- reading of non-standard branch ----- + for (Int_t iTrack = 0; iTrack < nTracksNonStd; iTrack++) { + AliVParticle *track = dynamic_cast ((*nonStdTracks)[iTrack]); + if (!track) + continue; + Double_t mom[3] = {track->Px(),track->Py(),track->Pz()}; p3.SetXYZ(mom[0],mom[1],mom[2]); pt = p3.Pt(); eta = p3.Eta(); - if (status == 0) continue; - if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue; - if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut + + // which cuts to apply if not AOD track (e.g. MC) ??? + AliAODTrack *trackAOD = dynamic_cast (track); + if (trackAOD) { + if (trackAOD->GetStatus() == 0) + continue; + if ((filterMask > 0) && !(trackAOD->TestFilterBit(filterMask))) + continue; + } + if ((eta > etaMax) || (eta < etaMin)) continue; // checking eta cut if (aodTrack == 0){ fRef->Delete(); // make sure to delete before placement new... @@ -308,18 +392,19 @@ Bool_t AliJetAODReader::FillMomentumArray() } new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag()); sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0; - cflag[aodTrack] = ( pt > ptMin ) ? 1: 0; + cflag[aodTrack] = ( pt > ptMin ) ? 1 : 0; aodTrack++; fRef->Add(track); + // printf("added non-standard track\n"); } - if(fDebug>0)printf("Used AOD tracks: %5d \n", aodTrack); + // set the signal flags fSignalFlag.Set(aodTrack,sflag); fCutFlag.Set(aodTrack,cflag); delete [] sflag; delete [] cflag; - + return kTRUE; } @@ -395,9 +480,11 @@ void AliJetAODReader::CreateTasks(TChain* tree) fDZ = fReaderHeader->GetDZ(); fTree = tree; - // Init EMCAL geometry and create UnitArray object - SetEMCALGeometry(); - // cout << "In create task" << endl; + // Init EMCAL geometry, if needed + if(fGeom == 0) + SetEMCALGeometry(); + else Info(" SetEMCALGeometry","was already done.. it's called just once !!"); + // Init parameters InitParameters(); InitUnitArray(); @@ -500,43 +587,14 @@ Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray) return kTRUE; } -//____________________________________________________________________________ -Bool_t AliJetAODReader::SetEMCALGeometry() -{ - // - // Set the EMCal Geometry - // - - if (!fTree->GetFile()) - return kFALSE; - TString geomFile(fTree->GetFile()->GetName()); - geomFile.ReplaceAll("AliESDs", "geometry"); - - // temporary workaround for PROOF bug #18505 - geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root"); - if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data()); - - // Define EMCAL geometry to be able to read ESDs - fGeom = AliJetDummyGeo::GetInstance(); - if (fGeom == 0) - fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL"); - - // To be setted to run some AliEMCALGeometry functions - TGeoManager::Import(geomFile); - fGeom->GetTransformationForSM(); - printf("\n EMCal Geometry set ! \n"); - - return kTRUE; - -} //____________________________________________________________________________ void AliJetAODReader::InitParameters() { // Initialise parameters fOpt = fReaderHeader->GetDetector(); - fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL + fNumUnits = fGeom->GetEMCGeometry()->GetNCells(); // Number of cells in EMCAL if(fDebug>1) printf("\n EMCal parameters initiated ! \n"); } @@ -573,31 +631,31 @@ void AliJetAODReader::InitUnitArray() // Define a grid of cell for the gaps between SM Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.; Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.; - fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0); + fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(0,phimin0,phimax0); fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015 fGrid0->SetGridType(0); fGrid0->SetMatrixIndexes(); fGrid0->SetIndexIJ(); n0 = fGrid0->GetNEntries(); - fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1); + fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(1,phimin1,phimax1); fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015 fGrid1->SetGridType(0); fGrid1->SetMatrixIndexes(); fGrid1->SetIndexIJ(); n1 = fGrid1->GetNEntries(); - fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2); + fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(2,phimin2,phimax2); fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015 fGrid2->SetGridType(0); fGrid2->SetMatrixIndexes(); fGrid2->SetIndexIJ(); n2 = fGrid2->GetNEntries(); - fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3); + fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(3,phimin3,phimax3); fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015 fGrid3->SetGridType(0); fGrid3->SetMatrixIndexes(); fGrid3->SetIndexIJ(); n3 = fGrid3->GetNEntries(); - fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4); + fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(4,phimin4,phimax4); fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015 fGrid4->SetGridType(0); fGrid4->SetMatrixIndexes();