X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliJetAODReader.cxx;h=1de092d38f7bc3b4d534f6be83eeba069d5eed1c;hb=c73e818954dddc41110e472137ac3e478d9f6782;hp=de57ba7afa42335dbed52f5fd0c66d30dbf39e47;hpb=f5c22dad61907c32ce04d1f1736147ecd353c244;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliJetAODReader.cxx b/JETAN/AliJetAODReader.cxx index de57ba7afa4..1de092d38f7 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 (AliEMCALGeoUtils) instead of dummy one (AliJetDummyGeo) +// moved geometry definition in AliJetReader +// marco.bregant@subatech.in2p3.fr //------------------------------------------------------------------------- @@ -31,16 +36,19 @@ #include #include #include +#include #include "AliJetAODReader.h" #include "AliJetAODReaderHeader.h" #include "AliAODEvent.h" #include "AliAODTrack.h" -#include "AliJetDummyGeo.h" +#include "AliAODMCParticle.h" +#include "AliEMCALGeoUtils.h" #include "AliJetAODFillUnitArrayTracks.h" #include "AliJetAODFillUnitArrayEMCalDigits.h" -#include "AliJetHadronCorrection.h" +#include "AliJetHadronCorrectionv1.h" #include "AliJetUnitArray.h" +#include "AliOADBContainer.h" ClassImp(AliJetAODReader) @@ -50,7 +58,6 @@ AliJetAODReader::AliJetAODReader(): fRef(new TRefArray), fDebug(0), fOpt(0), - fGeom(0), fHadCorr(0x0), fTpcGrid(0x0), fEmcalGrid(0x0), @@ -59,7 +66,6 @@ AliJetAODReader::AliJetAODReader(): fGrid2(0), fGrid3(0), fGrid4(0), - fPtCut(0), fApplyElectronCorrection(kFALSE), fApplyMIPCorrection(kTRUE), fApplyFractionHadronicCorrection(kFALSE), @@ -119,9 +125,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++; } } @@ -162,6 +166,98 @@ void AliJetAODReader::ConnectTree(TTree* tree, TObject* /*data*/) { } } + +Bool_t AliJetAODReader::AcceptAODMCParticle(AliAODMCParticle *mcP,Short_t flag){ + + // + // filter for charge and for charged and neutral, no detector + // response yet + // physical priamries already selected + + Int_t pdg = TMath::Abs(mcP->GetPdgCode()); + + // exclude neutrinos anyway + if((pdg == 12 || pdg == 14 || pdg == 16)) return kFALSE; + + if(flag==AliJetAODReaderHeader::kAllMC)return kTRUE; + if(flag==AliJetAODReaderHeader::kChargedMC){ + if(mcP->Charge()==0)return kFALSE; + return kTRUE; + } + + return kTRUE; + +} + + +Bool_t AliJetAODReader::FillMomentumArrayMC(){ + + // + // This routine fetches the MC particles from the AOD + // Depending on the flag all particles except neurinos are use + // or only charged particles + // + + TClonesArray *mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName()); + if(!mcArray){ + Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__); + return kFALSE; + } + + Int_t nMC = mcArray->GetEntriesFast(); + + // get number of tracks in event (for the loop) + if(fDebug>0)printf("AOD MC tracks: %5d \t", nMC); + + // temporary storage of signal and pt cut flag + Int_t* sflag = new Int_t[nMC]; + Int_t* cflag = new Int_t[nMC]; + + // get cuts set by user + Float_t ptMin = fReaderHeader->GetPtCut(); + Float_t etaMin = fReaderHeader->GetFiducialEtaMin(); + Float_t etaMax = fReaderHeader->GetFiducialEtaMax(); + Int_t mcTrack = 0; + Float_t pt, eta; + TVector3 p3; + + Short_t readerFlag = ((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC(); + + + for (Int_t it = 0; it < nMC; it++) { + AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it); + 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,readerFlag))continue; + pt = p3.Pt(); + + + + + if (mcTrack == 0){ + fRef->Delete(); // make sure to delete before placement new... + new(fRef) TRefArray(TProcessID::GetProcessWithUID(track)); + } + new ((*fMomentumArray)[mcTrack]) TLorentzVector(p3,p3.Mag()); + sflag[mcTrack] = 1; + cflag[mcTrack] = ( pt > ptMin ) ? 1: 0; + fRef->Add(track); + mcTrack++; + } + if(fDebug>0)printf("Used MC tracks: %5d \n", mcTrack); + // set the signal flags + fSignalFlag.Set(mcTrack,sflag); + fCutFlag.Set(mcTrack,cflag); + + delete [] sflag; + delete [] cflag; + + return kTRUE; +} + //____________________________________________________________________________ Bool_t AliJetAODReader::FillMomentumArray() @@ -170,56 +266,143 @@ Bool_t AliJetAODReader::FillMomentumArray() ClearArray(); fRef->Clear(); fDebug = fReaderHeader->GetDebug(); - + if (!fAOD) { return kFALSE; } - - // get number of tracks in event (for the loop) - Int_t nt = fAOD->GetNTracks(); - 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... + 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; + cflag[aodTrack] = ( pt > ptMin ) ? 1 : 0; aodTrack++; fRef->Add(track); + // printf("added non-standard track\n"); } - 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; } @@ -287,13 +470,19 @@ void AliJetAODReader::SetFractionHadronicCorrection(Double_t val) //____________________________________________________________________________ void AliJetAODReader::CreateTasks(TChain* tree) { + // + // For reader task initialization + // + fDebug = fReaderHeader->GetDebug(); 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(); @@ -322,7 +511,7 @@ void AliJetAODReader::CreateTasks(TChain* tree) fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection); fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection); fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection); - + // Add the task to global task fFillUnitArray->Add(fFillUAFromTracks); fFillUnitArray->Add(fFillUAFromEMCalDigits); fFillUAFromTracks->SetActive(kFALSE); @@ -348,6 +537,7 @@ Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray) // clear momentum array ClearArray(); + fRef->Clear(); fDebug = fReaderHeader->GetDebug(); fOpt = fReaderHeader->GetDetector(); @@ -362,6 +552,9 @@ Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray) fFillUAFromTracks->SetActive(kTRUE); fFillUAFromTracks->SetUnitArray(fUnitArray); fFillUAFromTracks->SetRefArray(fRefArray); + fFillUAFromTracks->SetReferences(fRef); + fFillUAFromTracks->SetSignalFlag(fSignalFlag); + fFillUAFromTracks->SetCutFlag(fCutFlag); fFillUAFromTracks->SetProcId(fProcId); // fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed fFillUAFromTracks->Exec("tpc"); @@ -369,6 +562,8 @@ Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray) fNumCandidate = fFillUAFromTracks->GetMult(); fNumCandidateCut = fFillUAFromTracks->GetMultCut(); } + fSignalFlag = fFillUAFromTracks->GetSignalFlag(); + fCutFlag = fFillUAFromTracks->GetCutFlag(); } // Digits only or Digits+TPC @@ -390,50 +585,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(); - // fHCorrection = 0; // For hadron correction - fHadCorr = 0; // For hadron correction - if(fEFlag==kFALSE){ - if(fOpt==0 || fOpt==1) - fECorrection = 0; // For electron correction - else fECorrection = 1; // For electron correction - } - 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"); } @@ -470,31 +629,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(); @@ -575,3 +734,4 @@ void AliJetAODReader::InitUnitArray() } // end grid type == 1 fArrayInitialised = 1; } +