// Jet ESD Reader
// ESD reader for jet analysis
// Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
-// Magali Estienne <magali.estienne@IReS.in2p3.fr>
+// Magali Estienne <magali.estienne@subatech.in2p3.fr>
//-------------------------------------------------------------------------
// --- Standard library ---
#include <TStopwatch.h>
#include <TLorentzVector.h>
#include <TVector3.h>
-#include <TTask.h>
+#include "TTask.h"
+#include "TTree.h"
+#include "TFile.h"
#include <TGeoManager.h>
#include <assert.h>
#include <TRefArray.h>
#include <TMath.h>
-#include <TChain.h>
-
+#include <TProcessID.h>
+#include <TRandom3.h>
// --- AliRoot header files ---
#include "AliJetESDReader.h"
#include "AliESD.h"
#include "AliESDtrack.h"
#include "AliJetDummyGeo.h"
-#include "AliJetFillUnitArrayTracks.h"
-#include "AliJetFillUnitArrayEMCalDigits.h"
+#include "AliJetESDFillUnitArrayTracks.h"
+#include "AliJetESDFillUnitArrayEMCalDigits.h"
+#include "AliJetHadronCorrectionv1.h"
#include "AliJetUnitArray.h"
+#include "AliAnalysisTask.h"
ClassImp(AliJetESDReader)
AliJetESDReader::AliJetESDReader():
AliJetReader(),
fGeom(0),
- fChain(0x0),
- fESD(0x0),
fHadCorr(0x0),
fTpcGrid(0x0),
fEmcalGrid(0x0),
fGrid3(0),
fGrid4(0),
fPtCut(0),
- fHCorrection(0),
+ fApplyElectronCorrection(kFALSE),
+ fApplyMIPCorrection(kTRUE),
+ fApplyFractionHadronicCorrection(kFALSE),
+ fFractionHadronicCorrection(0.3),
fNumUnits(0),
fDebug(0),
fMass(0),
fDZ(0),
fNeta(0),
fNphi(0),
- fArrayInitialised(0)
+ fRefArray(0x0),
+ fProcId(kFALSE)
{
// Constructor
}
{
// Destructor
delete fChain;
+ delete fTree;
delete fESD;
delete fTpcGrid;
delete fEmcalGrid;
if (strstr(name,pattern)){
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++;
int nMax = fChain->GetEntries();
printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);
-
+
// set number of events in header
if (fReaderHeader->GetLastEvent() == -1)
fReaderHeader->SetLastEvent(nMax);
Int_t nUsr = fReaderHeader->GetLastEvent();
fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
}
+
+}
+
+//__________________________________________________________
+void AliJetESDReader::SetApplyMIPCorrection(Bool_t val)
+{
+ //
+ // Set flag to apply MIP correction fApplyMIPCorrection
+ // - exclusive with fApplyFractionHadronicCorrection
+ //
+
+ fApplyMIPCorrection = val;
+ if(fApplyMIPCorrection == kTRUE)
+ {
+ SetApplyFractionHadronicCorrection(kFALSE);
+ printf("Enabling MIP Correction \n");
+ }
+ else
+ {
+ printf("Disabling MIP Correction \n");
+ }
+}
+
+//__________________________________________________________
+void AliJetESDReader::SetApplyFractionHadronicCorrection(Bool_t val)
+{
+ //
+ // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
+ // - exclusive with fApplyMIPCorrection
+ //
+
+ fApplyFractionHadronicCorrection = val;
+ if(fApplyFractionHadronicCorrection == kTRUE)
+ {
+ SetApplyMIPCorrection(kFALSE);
+ printf("Enabling Fraction Hadronic Correction \n");
+ }
+ else
+ {
+ printf("Disabling Fraction Hadronic Correction \n");
+ }
+}
+
+//__________________________________________________________
+void AliJetESDReader::SetFractionHadronicCorrection(Double_t val)
+{
+ //
+ // Set value to fFractionHadronicCorrection (default is 0.3)
+ // apply EMC hadronic correction fApplyFractionHadronicCorrection
+ // - exclusive with fApplyMIPCorrection
+ //
+
+ fFractionHadronicCorrection = val;
+ if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
+ {
+ SetApplyFractionHadronicCorrection(kTRUE);
+ printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
+ }
+ else
+ {
+ SetApplyFractionHadronicCorrection(kFALSE);
+ }
}
+
//____________________________________________________________________________
-void AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) {
+void AliJetESDReader::SetInputEvent(const TObject* esd, const TObject* /*aod*/, const TObject* /*mc*/) {
// Connect the tree
fESD = (AliESDEvent*) esd;
}
//____________________________________________________________________________
-Bool_t AliJetESDReader::FillMomentumArray(Int_t /*event*/)
+Bool_t AliJetESDReader::FillMomentumArray()
{
// Fill momentum array
// 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];
Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
- //loop over tracks
+ //loop over tracks in ESD
for (Int_t it = 0; it < nt; it++) {
AliESDtrack *track = fESD->GetTrack(it);
UInt_t status = track->GetStatus();
if (pt > ptMin) cflag[goodTrack]=1; // pt cut
goodTrack++;
}
+
// set the signal flags
fSignalFlag.Set(goodTrack,sflag);
fCutFlag.Set(goodTrack,cflag);
+
+ delete[] sflag;
+ delete[] cflag;
+
return kTRUE;
}
//____________________________________________________________________________
-void AliJetESDReader::CreateTasks()
+void AliJetESDReader::CreateTasks(TChain* tree)
{
+ //
+ // For reader task initialization
+ //
+
fDebug = fReaderHeader->GetDebug();
fDZ = fReaderHeader->GetDZ();
+ fTree = tree;
- // Init EMCAL geometry and create UnitArray object
+ // 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");
- fFillUAFromTracks = new AliJetFillUnitArrayTracks();
+ // Create a task for to fill the charged particle information
+ fFillUAFromTracks = new AliJetESDFillUnitArrayTracks();
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->SetApplyMIPCorrection(fApplyMIPCorrection);
fFillUAFromTracks->SetHadCorrector(fHadCorr);
- fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits();
+ // Create a task for to fill the neutral particle information
+ fFillUAFromEMCalDigits = new AliJetESDFillUnitArrayEMCalDigits();
fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
fFillUAFromEMCalDigits->SetGeom(fGeom);
fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
- // fFillUnitArray->Add(fFillUAFromTracks);
+ fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
+ fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
+ fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
+ // Add the task to global task
+ fFillUnitArray->Add(fFillUAFromTracks);
fFillUnitArray->Add(fFillUAFromEMCalDigits);
fFillUAFromTracks->SetActive(kFALSE);
fFillUAFromEMCalDigits->SetActive(kFALSE);
}
//____________________________________________________________________________
-//void AliJetESDReader::ExecTasks(Int_t event)
-Bool_t AliJetESDReader::ExecTasks(Int_t /*event*/)
+Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray)
{
- // clear momentum array
- Int_t nEntRef = fRefArray->GetEntries();
-
- for(Int_t i=0; i<nEntRef; i++)
- {
- ((AliJetUnitArray*)fRefArray->At(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);
- }
+ //
+ // Reader task execussion
+ //
+
+ fProcId = procid;
+ fRefArray = refArray;
+ // clear momentum array
ClearArray();
fDebug = fReaderHeader->GetDebug();
fOpt = fReaderHeader->GetDetector();
- // InitParameters();
if(!fESD) {
return kFALSE;
}
- /*
- // 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");
+ fFillUAFromTracks->SetProcId(fProcId);
+ // fFillUAFromTracks->ExecuteTask("tpc"); // Temporarily changed
+ fFillUAFromTracks->Exec("tpc"); // Temporarily added
if(fOpt==1){
fNumCandidate = fFillUAFromTracks->GetMult();
fNumCandidateCut = fFillUAFromTracks->GetMultCut();
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 changed !!!
+ fFillUAFromEMCalDigits->Exec("digits"); // Temporarily added
fNumCandidate = fFillUAFromEMCalDigits->GetMult();
fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
}
- // fFillUnitArray->ExecuteTask(); // => Temporarily commented
+ // fFillUnitArray->ExecuteTask(); // Temporarily commented
return kTRUE;
}
//____________________________________________________________________________
-void AliJetESDReader::SetEMCALGeometry()
+Bool_t 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");
+ //
+ // Set the EMCal Geometry
+ //
- // To be setted to run some AliEMCALGeometry functions
- TGeoManager::Import("geometry.root");
- // fGeom->GetTransformationForSM();
- printf("\n EMCal Geometry set ! \n");
+ 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
- fHCorrection = 0; // For hadron correction
- fHadCorr = 0; // For hadron correction
- fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL
- if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
+ // 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()
{
- //Initialises unit arrays
+ //
+ // Create and Initialise unit arrays
+ //
+
Int_t nElements = fTpcGrid->GetNEntries();
Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
if(fArrayInitialised) fUnitArray->Delete();
// 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);
+ new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
}
}
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;
}
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);
+ new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
}
else {
if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
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);
+ 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 && nBin<fNumUnits+nElements+nGaps){
if(nBin<fNumUnits+nElements+n0)
{
- Float_t phi = eta = 0.;
- fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
+ Float_t phi0 = eta = 0.;
+ fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi0,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);
+ new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
}
else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
{
- Float_t phi = eta = 0.;
- fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
+ Float_t phi0 = eta = 0.;
+ fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi0,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);
+ 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 && nBin<fNumUnits+nElements+n0+n1+n2)
{
- Float_t phi = eta = 0.;
- fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
+ Float_t phi0 = eta = 0.;
+ fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi0,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);
+ 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 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
{
- Float_t phi = eta = 0.;
- fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
+ Float_t phi0 = eta = 0.;
+ fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi0,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);
+ 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 && nBin<fNumUnits+nElements+nGaps)
{
- Float_t phi = eta = 0.;
- fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
+ Float_t phi0 = eta = 0.;
+ fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi0,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);
+ new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
}
}
} // end if(fDZ)
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff