X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliJetESDReader.cxx;h=1ac583d8f6ac3245faf2c53d1ec574a3b86e3561;hb=e89c03c6eef54e1640506dc66a117fbcc7f9d369;hp=318b37cf3baf8e4bfcfc7a3d23f663fd654e9af2;hpb=69b186410c2bd37df6d4d6460607f2e7b7762235;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliJetESDReader.cxx b/JETAN/AliJetESDReader.cxx index 318b37cf3ba..1ac583d8f6a 100755 --- a/JETAN/AliJetESDReader.cxx +++ b/JETAN/AliJetESDReader.cxx @@ -20,48 +20,65 @@ // Magali Estienne //------------------------------------------------------------------------- - +// --- Standard library --- #include + +// --- ROOT system --- #include +#include #include #include +#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 "AliEMCALGeometry.h" #include "AliJetDummyGeo.h" #include "AliJetFillUnitArrayTracks.h" +#include "AliJetFillUnitArrayEMCalDigits.h" #include "AliJetUnitArray.h" ClassImp(AliJetESDReader) AliJetESDReader::AliJetESDReader(): - AliJetReader(), - fGeom(0), - fChain(0x0), - fESD(0x0), - fHadCorr(0x0), - fTpcGrid(0x0), - fEmcalGrid(0x0), - fPtCut(0), - fHCorrection(0), - fNumUnits(0), - fDebug(0), - fNIn(0), - fOpt(0), - fNeta(0), - fNphi(0), - fArrayInitialised(0) + AliJetReader(), + fGeom(0), + fChain(0x0), + fESD(0x0), + fHadCorr(0x0), + fTpcGrid(0x0), + fEmcalGrid(0x0), + fGrid0(0), + fGrid1(0), + fGrid2(0), + fGrid3(0), + fGrid4(0), + fPtCut(0), + fHCorrection(0), + fNumUnits(0), + fDebug(0), + fMass(0), + fSign(0), + fNIn(0), + fOpt(0), + fDZ(0), + fNeta(0), + fNphi(0), + fArrayInitialised(0) { - // Constructor - SetEMCALGeometry(); + // Constructor } //____________________________________________________________________________ - AliJetESDReader::~AliJetESDReader() { // Destructor @@ -69,10 +86,17 @@ AliJetESDReader::~AliJetESDReader() 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 @@ -82,8 +106,7 @@ void AliJetESDReader::OpenInputFiles() const char* dirName=fReaderHeader->GetDirectory(); const char* pattern=fReaderHeader->GetPattern(); -// // Add files matching patters to the chain - + // Add files matching patters to the chain void *dir = gSystem->OpenDirectory(dirName); const char *name = 0x0; int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd(); @@ -92,22 +115,23 @@ void AliJetESDReader::OpenInputFiles() if (a>=nesd) continue; if (strstr(name,pattern)){ - char path[256]; - sprintf(path,"%s/%s/AliESDs.root",dirName,name); - fChain->AddFile(path); - a++; + printf("Adding %s\n",name); + char path[256]; + // sprintf(path,"%s/%s/AliESDs.root",dirName,name); + sprintf(path,"%s/%s/AliESDs.root",dirName,name); + fChain->AddFile(path); + a++; } } gSystem->FreeDirectory(dir); + fESD = new AliESDEvent(); + fESD->ReadFromTree(fChain); - fESD = 0; - fChain->SetBranchAddress("ESD", &fESD); - 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) @@ -118,25 +142,14 @@ void AliJetESDReader::OpenInputFiles() } } -void AliJetESDReader::ConnectTree(TTree* tree) { +//____________________________________________________________________________ +void AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) { // Connect the tree - fChain = (TChain*) tree; - - fChain->SetBranchAddress("ESD", &fESD); - Int_t nMax = fChain->GetEntries(); - printf("\nTotal number of events in chain= %5d", nMax); - // set number of events in header - if (fReaderHeader->GetLastEvent() == -1) - fReaderHeader->SetLastEvent(nMax); - else { - Int_t nUsr = fReaderHeader->GetLastEvent(); - fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr)); - } + fESD = (AliESDEvent*) esd; } //____________________________________________________________________________ - -Bool_t AliJetESDReader::FillMomentumArray(Int_t event) +Bool_t AliJetESDReader::FillMomentumArray() { // Fill momentum array @@ -148,15 +161,16 @@ Bool_t AliJetESDReader::FillMomentumArray(Int_t event) // clear momentum array ClearArray(); fDebug = fReaderHeader->GetDebug(); - // get event from chain - fChain->GetTree()->GetEntry(event); + 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]; @@ -195,49 +209,140 @@ Bool_t AliJetESDReader::FillMomentumArray(Int_t event) fSignalFlag.Set(goodTrack,sflag); fCutFlag.Set(goodTrack,cflag); -// -// - if (fTpcGrid || fEmcalGrid) { - SetEMCALGeometry(); - InitParameters(); - AliJetFillUnitArrayTracks *fillUAFromTracks = new AliJetFillUnitArrayTracks(); - fillUAFromTracks->SetReaderHeader(fReaderHeader); - fillUAFromTracks->SetMomentumArray(fMomentumArray); - fillUAFromTracks->SetTPCGrid(fTpcGrid); - fillUAFromTracks->SetEMCalGrid(fEmcalGrid); - fillUAFromTracks->SetHadCorrection(fHCorrection); - fillUAFromTracks->SetHadCorrector(fHadCorr); - fNeta = fillUAFromTracks->GetNeta(); - fNphi = fillUAFromTracks->GetNphi(); - fillUAFromTracks->SetActive(kFALSE); - // TPC only or Digits+TPC or Clusters+TPC - if(fOpt%2==!0 && fOpt!=0) { - fillUAFromTracks->SetActive(kTRUE); - fillUAFromTracks->SetUnitArray(fUnitArray); - fillUAFromTracks->ExecuteTask("tpc"); - } + delete[] sflag; + delete[] cflag; + + return kTRUE; +} + +//____________________________________________________________________________ +void AliJetESDReader::CreateTasks() +{ + fDebug = fReaderHeader->GetDebug(); + fDZ = fReaderHeader->GetDZ(); + + // Init EMCAL geometry and create UnitArray object + SetEMCALGeometry(); + InitParameters(); + InitUnitArray(); + + fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder"); + fFillUAFromTracks = new AliJetFillUnitArrayTracks(); + fFillUAFromTracks->SetReaderHeader(fReaderHeader); + fFillUAFromTracks->SetGeom(fGeom); + fFillUAFromTracks->SetTPCGrid(fTpcGrid); + fFillUAFromTracks->SetEMCalGrid(fEmcalGrid); + + if(fDZ) + { + fFillUAFromTracks->SetGrid0(fGrid0); + fFillUAFromTracks->SetGrid1(fGrid1); + fFillUAFromTracks->SetGrid2(fGrid2); + fFillUAFromTracks->SetGrid3(fGrid3); + fFillUAFromTracks->SetGrid4(fGrid4); + } + fFillUAFromTracks->SetHadCorrection(fHCorrection); + fFillUAFromTracks->SetHadCorrector(fHadCorr); + fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits(); + fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader); + fFillUAFromEMCalDigits->SetGeom(fGeom); + fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid); + fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid); + // 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; + +} + +//____________________________________________________________________________ +//void AliJetESDReader::ExecTasks(Int_t event) +Bool_t AliJetESDReader::ExecTasks(Int_t /*event*/) +{ + // clear momentum array + Int_t nEntRef = fRefArray->GetEntries(); + + for(Int_t i=0; iAt(i))->SetUnitTrackID(0); + ((AliJetUnitArray*)fRefArray->At(i))->SetUnitEnergy(0.); + ((AliJetUnitArray*)fRefArray->At(i))->SetUnitCutFlag(kPtSmaller); + ((AliJetUnitArray*)fRefArray->At(i))->SetUnitDetectorFlag(kTpc); + ((AliJetUnitArray*)fRefArray->At(i))->SetUnitFlag(kOutJet); + } + + ClearArray(); + + fDebug = fReaderHeader->GetDebug(); + fOpt = fReaderHeader->GetDetector(); + // InitParameters(); + + if(!fESD) { + return kFALSE; + } - delete fillUAFromTracks; + /* + // get event from chain + // For TSelectors + // fChain->GetTree()->GetEntry(event); + // For interactive process + // fChain->GetEntry(event); + fChain->GetEvent(event); + */ + + // 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->ExecuteTask("tpc"); // => Temporarily changed !!! + // Incompatibility with Andrei's analysis framework + fFillUAFromTracks->Exec("tpc"); + 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->SetInitMult(fFillUAFromTracks->GetMult()); + fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut()); + fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily changed !!! + fNumCandidate = fFillUAFromEMCalDigits->GetMult(); + fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut(); + } + + // fFillUnitArray->ExecuteTask(); // => Temporarily commented + return kTRUE; } - +//____________________________________________________________________________ void AliJetESDReader::SetEMCALGeometry() { // Define EMCAL geometry to be able to read ESDs fGeom = AliJetDummyGeo::GetInstance(); if (fGeom == 0) fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL"); - + // To be setted to run some AliEMCALGeometry functions TGeoManager::Import("geometry.root"); - fGeom->GetTransformationForSM(); + // fGeom->GetTransformationForSM(); printf("\n EMCal Geometry set ! \n"); } - + + +//____________________________________________________________________________ void AliJetESDReader::InitParameters() { // Initialise parameters @@ -247,16 +352,153 @@ void AliJetESDReader::InitParameters() if(fDebug>1) printf("\n EMCal parameters initiated ! \n"); } +//____________________________________________________________________________ void AliJetESDReader::InitUnitArray() { //Initialises unit arrays - Int_t nElements = fTpcGrid->GetNEntries(); - if(fArrayInitialised) delete [] fUnitArray; - if(fTpcGrid->GetGridType()==0) - fUnitArray = new AliJetUnitArray[nElements]; - if(fTpcGrid->GetGridType()==1) - fUnitArray = new AliJetUnitArray[fNumUnits+nElements]; - fArrayInitialised = 1; + 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,eta,phi); + 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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + } + } + + if(fTpcGrid->GetGridType()==1) + { + Int_t nGaps = 0; + Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0; + + 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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + } + else { + if(fDZ) { + if(nBin>=fNumUnits+nElements && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta); + Deta = fGrid0->GetDeta(); + Dphi = fGrid0->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta); + Deta = fGrid1->GetDeta(); + Dphi = fGrid1->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0+n1 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta); + Deta = fGrid2->GetDeta(); + Dphi = fGrid2->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta); + Deta = fGrid3->GetDeta(); + Dphi = fGrid3->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta); + Deta = fGrid4->GetDeta(); + Dphi = fGrid4->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + } + } + } // end if(fDZ) + } // end else 2 + } // end else 1 + } // end loop on nBin + } // end grid type == 1 + fArrayInitialised = 1; } +