X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=inline;f=JETAN%2FAliJetESDReader.cxx;h=6675a5d5b48e1c9ac9dcf74571ee4ce3032c56ae;hb=b682a4ecf8d555d2a2c68f449157e928e172d8bf;hp=956bb063ae765af6737edb14d68ae4f6b481497a;hpb=99e5fe42da88af9a55105ef427fcea810c10b194;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliJetESDReader.cxx b/JETAN/AliJetESDReader.cxx index 956bb063ae7..6675a5d5b48 100755 --- a/JETAN/AliJetESDReader.cxx +++ b/JETAN/AliJetESDReader.cxx @@ -12,75 +12,134 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - -// + +//------------------------------------------------------------------------- // Jet ESD Reader // ESD reader for jet analysis -// Author: Mercedes Lopez Noriega -// mercedes.lopez.noriega@cern.ch -// - +// Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch) +// Magali Estienne +//------------------------------------------------------------------------- +// --- Standard library --- #include + +// --- ROOT system --- #include +#include #include +#include +#include "TTask.h" +#include "TTree.h" +#include "TFile.h" +#include +#include +#include +#include +#include +#include + +// --- AliRoot header files --- #include "AliJetESDReader.h" #include "AliJetESDReaderHeader.h" +#include "AliESDEvent.h" #include "AliESD.h" #include "AliESDtrack.h" +#include "AliJetDummyGeo.h" +#include "AliJetFillUnitArrayTracks.h" +#include "AliJetFillUnitArrayEMCalDigits.h" +#include "AliJetUnitArray.h" +#include "AliAnalysisTask.h" ClassImp(AliJetESDReader) -AliJetESDReader::AliJetESDReader() +AliJetESDReader::AliJetESDReader(): + AliJetReader(), + fGeom(0), + fChain(0x0), + fTree(0x0), + fESD(0x0), + fHadCorr(0x0), + fTpcGrid(0x0), + fEmcalGrid(0x0), + fGrid0(0), + fGrid1(0), + fGrid2(0), + fGrid3(0), + fGrid4(0), + fPtCut(0), + fHCorrection(0), + fECorrection(0), + fEFlag(kFALSE), + fNumUnits(0), + fDebug(0), + fMass(0), + fSign(0), + fNIn(0), + fOpt(0), + fDZ(0), + fNeta(0), + fNphi(0), + fArrayInitialised(0), + fRefArray(0x0), + fProcId(kFALSE) { - // Constructor - - fReaderHeader = 0x0; - fMass = 0; - fSign = 0; + // Constructor } //____________________________________________________________________________ - AliJetESDReader::~AliJetESDReader() { // Destructor - delete fReaderHeader; + delete fChain; + delete fTree; + delete fESD; + delete fTpcGrid; + delete fEmcalGrid; + if(fDZ) + { + delete fGrid0; + delete fGrid1; + delete fGrid2; + delete fGrid3; + delete fGrid4; + } } //____________________________________________________________________________ - void AliJetESDReader::OpenInputFiles() - { // chain for the ESDs fChain = new TChain("esdTree"); - fChainMC = new TChain("mcStackTree"); // get directory and pattern name from the header - const char* dirName=fReaderHeader->GetDirectory(); - const char* pattern=fReaderHeader->GetPattern(); + const char* dirName=fReaderHeader->GetDirectory(); + const char* pattern=fReaderHeader->GetPattern(); - // Add files matching patters to the chain - void *dir = gSystem->OpenDirectory(dirName); - const char *name = 0x0; - while ((name = gSystem->GetDirEntry(dir))){ - if (strstr(name,pattern)){ - printf("Adding %s\n",name); - char path[256]; - sprintf(path,"%s/%s",dirName,name); - fChain->AddFile(path,-1,"esdTree"); - fChainMC->AddFile(path,-1,"mcStackTree"); - } - } - - gSystem ->FreeDirectory(dir); - fChain ->SetBranchAddress("ESD", &fESD); - fChainMC->SetBranchAddress("Header", &fAliHeader); - fChainMC->SetBranchAddress("Stack", &fArrayMC); + // Add files matching patters to the chain + void *dir = gSystem->OpenDirectory(dirName); + const char *name = 0x0; + int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd(); + int a = 0; + while ((name = gSystem->GetDirEntry(dir))){ + if (a>=nesd) continue; + + if (strstr(name,pattern)){ + printf("Adding %s\n",name); + char path[256]; + sprintf(path,"%s/%s/AliESDs.root",dirName,name); + fChain->AddFile(path); + a++; + } + } + + gSystem->FreeDirectory(dir); + + fESD = new AliESDEvent(); + fESD->ReadFromTree(fChain); int nMax = fChain->GetEntries(); - printf("\nTotal number of events in chain= %d",nMax); + + printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax); // set number of events in header if (fReaderHeader->GetLastEvent() == -1) @@ -89,48 +148,387 @@ void AliJetESDReader::OpenInputFiles() Int_t nUsr = fReaderHeader->GetLastEvent(); fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr)); } + } //____________________________________________________________________________ +void AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) { + // Connect the tree + fESD = (AliESDEvent*) esd; +} -void AliJetESDReader::FillMomentumArray(Int_t event) +//____________________________________________________________________________ +Bool_t AliJetESDReader::FillMomentumArray() { -// Fills the momentum array + // Fill momentum array + Int_t goodTrack = 0; Int_t nt = 0; - Float_t pt; - Double_t p[3]; // track 3 momentum vector - + Float_t pt, eta; + TVector3 p3; + // clear momentum array ClearArray(); - // get event from chain - fChain->GetEntry(event); - fChainMC->GetEntry(event); + fDebug = fReaderHeader->GetDebug(); + fOpt = fReaderHeader->GetDetector(); + + if (!fESD) { + return kFALSE; + } + // get number of tracks in event (for the loop) nt = fESD->GetNumberOfTracks(); + printf("Fill Momentum Array %5d \n", 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 = ((AliJetESDReaderHeader*) fReaderHeader)->GetPtCut(); - - //loop over tracks + Float_t ptMin = fReaderHeader->GetPtCut(); + Float_t etaMin = fReaderHeader->GetFiducialEtaMin(); + Float_t etaMax = fReaderHeader->GetFiducialEtaMax(); + + //loop over tracks in ESD for (Int_t it = 0; it < nt; it++) { AliESDtrack *track = fESD->GetTrack(it); UInt_t status = track->GetStatus(); - if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check + + Double_t mom[3]; + track->GetPxPyPz(mom); + p3.SetXYZ(mom[0],mom[1],mom[2]); + pt = p3.Pt(); + if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check + if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check + if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() + && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check + if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly() + && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check + eta = p3.Eta(); + if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut + + new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag()); + sflag[goodTrack]=0; + if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; + cflag[goodTrack]=0; + if (pt > ptMin) cflag[goodTrack]=1; // pt cut + goodTrack++; + } - if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() && TMath::Abs(track->GetLabel()) > 10000) continue; - - track->GetPxPyPz(p); - pt = TMath::Sqrt(p[0] * p[0] + p[1] * p[1]); // pt of the track + // set the signal flags + fSignalFlag.Set(goodTrack,sflag); + fCutFlag.Set(goodTrack,cflag); - if (pt < ptMin) continue; //check cuts + delete[] sflag; + delete[] cflag; - new ((*fMomentumArray)[goodTrack]) - TLorentzVector(p[0], p[1], p[2], - TMath::Sqrt(pt * pt +p[2] * p[2])); - goodTrack++; + return kTRUE; + +} + +//____________________________________________________________________________ +void AliJetESDReader::CreateTasks(TChain* tree) +{ + // + // For reader task initialization + // + + fDebug = fReaderHeader->GetDebug(); + fDZ = fReaderHeader->GetDZ(); + fTree = tree; + + // Init EMCAL geometry + SetEMCALGeometry(); + // Init parameters + InitParameters(); + // Create and init unit array + InitUnitArray(); + + // Create global reader task for analysis + fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder"); + // Create a task for to fill the charged particle information + fFillUAFromTracks = new AliJetFillUnitArrayTracks(); + fFillUAFromTracks->SetReaderHeader(fReaderHeader); + fFillUAFromTracks->SetGeom(fGeom); + fFillUAFromTracks->SetTPCGrid(fTpcGrid); + fFillUAFromTracks->SetEMCalGrid(fEmcalGrid); + if(fDZ) + { // Calo dead zones inclusion + fFillUAFromTracks->SetGrid0(fGrid0); + fFillUAFromTracks->SetGrid1(fGrid1); + fFillUAFromTracks->SetGrid2(fGrid2); + fFillUAFromTracks->SetGrid3(fGrid3); + fFillUAFromTracks->SetGrid4(fGrid4); + } + fFillUAFromTracks->SetHadCorrection(fHCorrection); + fFillUAFromTracks->SetHadCorrector(fHadCorr); + // Create a task for to fill the neutral particle information + fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits(); + fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader); + fFillUAFromEMCalDigits->SetGeom(fGeom); + fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid); + fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid); + fFillUAFromEMCalDigits->SetEleCorrection(fECorrection); + // Add the task to global task + fFillUnitArray->Add(fFillUAFromTracks); + fFillUnitArray->Add(fFillUAFromEMCalDigits); + fFillUAFromTracks->SetActive(kFALSE); + fFillUAFromEMCalDigits->SetActive(kFALSE); + + cout << "Tasks instantiated at that stage ! " << endl; + cout << "You can loop over events now ! " << endl; + +} + +//____________________________________________________________________________ +Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray) +{ + // + // Reader task execussion + // + + fProcId = procid; + fRefArray = refArray; + vector vtmp(3); + + // clear momentum array + ClearArray(); + + fDebug = fReaderHeader->GetDebug(); + fOpt = fReaderHeader->GetDetector(); + + if(!fESD) { + return kFALSE; + } + + // TPC only or Digits+TPC or Clusters+TPC + if(fOpt%2==!0 && fOpt!=0){ + fFillUAFromTracks->SetESD(fESD); + fFillUAFromTracks->SetActive(kTRUE); + fFillUAFromTracks->SetUnitArray(fUnitArray); + fFillUAFromTracks->SetRefArray(fRefArray); + fFillUAFromTracks->SetProcId(fProcId); + // fFillUAFromTracks->ExecuteTask("tpc"); // Temporarily changed + fFillUAFromTracks->Exec("tpc"); // Temporarily added + if(fOpt==1){ + fNumCandidate = fFillUAFromTracks->GetMult(); + fNumCandidateCut = fFillUAFromTracks->GetMultCut(); + } + } + + // Digits only or Digits+TPC + if(fOpt>=2 && fOpt<=3){ + fFillUAFromEMCalDigits->SetESD(fESD); + fFillUAFromEMCalDigits->SetActive(kTRUE); + fFillUAFromEMCalDigits->SetUnitArray(fUnitArray); + fFillUAFromEMCalDigits->SetRefArray(fRefArray); + fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId()); + fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult()); + fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut()); + fFillUAFromEMCalDigits->Exec("digits"); // Temporarily added + fNumCandidate = fFillUAFromEMCalDigits->GetMult(); + fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut(); + } + + // fFillUnitArray->ExecuteTask(); // Temporarily commented + + return kTRUE; +} + +//____________________________________________________________________________ +Bool_t AliJetESDReader::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 AliJetESDReader::InitParameters() +{ + // Initialise parameters + fOpt = fReaderHeader->GetDetector(); + 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 + if(fDebug>1) printf("\n EMCal parameters initiated ! \n"); +} + +//____________________________________________________________________________ +void AliJetESDReader::InitUnitArray() +{ + // + // Create and Initialise unit arrays + // + + Int_t nElements = fTpcGrid->GetNEntries(); + Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.; + if(fArrayInitialised) fUnitArray->Delete(); + + if(fTpcGrid->GetGridType()==0) + { // Fill the following quantities : + // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi, + // detector flag, in/out jet, pt cut, mass, cluster ID) + for(Int_t nBin = 1; nBin < nElements+1; nBin++) + { + fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta); + phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi); + Deta = fTpcGrid->GetDeta(); + Dphi = fTpcGrid->GetDphi(); + new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + } + + if(fTpcGrid->GetGridType()==1) + { + Int_t nGaps = 0; + Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0; - printf("\nIn event %d, number of good tracks %d \n", event, goodTrack); + if(fDZ) + { + // 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); + 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); + 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); + 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); + 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); + fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015 + fGrid4->SetGridType(0); + fGrid4->SetMatrixIndexes(); + fGrid4->SetIndexIJ(); + n4 = fGrid4->GetNEntries(); + + if(fDebug>1) + { + cout << "n0 cells: " << n0 << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl; + cout << "n1 cells: " << n1 << "phimin1: " << phimin1 << ", phimax1: " << phimax1 << endl; + cout << "n2 cells: " << n2 << "phimin2: " << phimin2 << ", phimax2: " << phimax2 << endl; + cout << "n3 cells: " << n3 << "phimin3: " << phimin3 << ", phimax3: " << phimax3 << endl; + cout << "n4 cells: " << n4 << "phimin4: " << phimin4 << ", phimax4: " << phimax4 << endl; + } + + nGaps = n0+n1+n2+n3+n4; + + } + + for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++) + { + if(nBinEtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry + // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid + phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi); + Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values + Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else { + if(nBin>=fNumUnits && nBinGetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta); + phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi); + Deta = fTpcGrid->GetDeta(); + Dphi = fTpcGrid->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else { + if(fDZ) { + if(nBin>=fNumUnits+nElements && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi0,eta); + Deta = fGrid0->GetDeta(); + Dphi = fGrid0->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi0,eta); + Deta = fGrid1->GetDeta(); + Dphi = fGrid1->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0+n1 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi0,eta); + Deta = fGrid2->GetDeta(); + Dphi = fGrid2->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi0,eta); + Deta = fGrid3->GetDeta(); + Dphi = fGrid3->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi0,eta); + Deta = fGrid4->GetDeta(); + Dphi = fGrid4->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + } + } // end if(fDZ) + } // end else 2 + } // end else 1 + } // end loop on nBin + } // end grid type == 1 + fArrayInitialised = 1; } +