From e36a3f229e62a74f64da011aed64fe003535eaa8 Mon Sep 17 00:00:00 2001 From: morsch Date: Fri, 2 Oct 2009 08:27:57 +0000 Subject: [PATCH] Improved hadron correction (Magali, Mark) Track references in rec. Jet for UA1v2 (Magali) --- JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx | 135 ++++++++++++++-- JETAN/AliJetAODFillUnitArrayEMCalDigits.h | 4 +- JETAN/AliJetAODFillUnitArrayTracks.cxx | 112 ++++++++++--- JETAN/AliJetAODFillUnitArrayTracks.h | 2 +- JETAN/AliJetAODReader.cxx | 17 +- JETAN/AliJetESDFillUnitArrayTracks.cxx | 34 +++- JETAN/AliJetESDReader.cxx | 5 +- JETAN/AliJetFillUnitArray.cxx | 3 + JETAN/AliJetFillUnitArray.h | 11 ++ JETAN/AliJetHadronCorrectionv0.cxx | 8 +- JETAN/AliJetHadronCorrectionv0.h | 3 +- JETAN/AliJetHadronCorrectionv1.cxx | 170 ++++++++++++-------- JETAN/AliJetHadronCorrectionv1.h | 43 ++--- JETAN/AliUA1JetFinderV2.cxx | 25 +++ 14 files changed, 423 insertions(+), 149 deletions(-) diff --git a/JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx b/JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx index 2e9429d913b..0b4151d54ba 100755 --- a/JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx +++ b/JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx @@ -138,18 +138,13 @@ void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/) fCluster = fReaderHeader->GetCluster(); // Init parameters - InitParameters(); + // InitParameters(); -//(not used ?) Int_t nDigitTot = 0; Int_t goodDigit = 0; -//(not used ?) Int_t beg = 0; -//(not used ?) Int_t end = 0; Int_t index = 0; -//(not used ?) Int_t count = 0; - -//(not used ?) Int_t nRefEnt = fRefArray->GetEntries(); if(!fCluster) { // Keep all digit information + // Loop over all cell information //------------------------------------------------------------------ AliAODCaloCells &cells= *(fAOD->GetEMCALCells()); @@ -168,7 +163,7 @@ void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/) phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD); - Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); + // Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID); @@ -189,6 +184,46 @@ void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/) ok = kTRUE; } + // Hadronic Correction + double correction=0; + if (fApplyFractionHadronicCorrection) { + + TArrayS* matched = new TArrayS(); + GetTracksPointingToCell(matched, etaD, phiD,0.02); + + double ptot = 0; + for(int itrack = 0; itrack < matched->GetSize(); itrack++) + { + const short index = matched->At(itrack); + if (index>-1) + { + AliAODTrack *track = fAOD->GetTrack(index); + ptot += track->P(); + //printf("aod-track p=%f \n",track->P() ); + } + } + + correction = ptot * fFractionHadronicCorrection; + + if (ptot>0 && fDebug>1 ) { + printf("SUM of track momentum for this cell: p=%f \n", ptot); + printf("fractional correction=%f \n", fFractionHadronicCorrection); + printf("TOTAL correction=%f \n", correction ); + } + + delete matched; + + }//end hadronic correction + + double enerCorr = digitEn; + if (correction >= enerCorr) enerCorr = 0; + else enerCorr -= correction; + if (correction>0 && fDebug>1) { + printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr); + } + + Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); + // Detector flag if(unitEnergy>0.) uArray->SetUnitDetectorFlag(kAll); @@ -224,26 +259,22 @@ void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/) // Get CaloCells AliAODCaloCells &cells= *(fAOD->GetEMCALCells()); -//(not used ?) Int_t ncell = cells.GetNumberOfCells() ; -//(not used ?) Int_t type = cells.GetType(); for(Int_t j = beg; j < nclus; j++) { // loop over clusters // Retrieve cluster from aod AliAODCaloCluster *fClus = (AliAODCaloCluster *) caloClusters->At(j) ; // Get the cluster info -//(not used ?) Float_t energy = fClus->E() ; fClus->GetPosition(pos) ; TVector3 vpos(pos[0],pos[1],pos[2]) ; Int_t digMult = fClus->GetNCells() ; UShort_t *digID = fClus->GetCellsAbsId() ; -//(not used ?) Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ; Int_t trackIndex = -1; if(fClus->GetNTracksMatched()!=0) trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID(); -//(not used ?) Int_t cLabel = fClus->GetLabel(0); + // Do double-counted electron correction if (fApplyElectronCorrection != 0 && trackIndex !=-1 ) { @@ -254,7 +285,6 @@ void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/) // Get CaloCells of cluster and fill the unitArray for(Int_t i = 0; i < digMult ; i++) { Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ; -//(not used ?) Double_t digitAmpFrac = digAmpFrac[i]; Float_t digitAmp = cells.GetCellAmplitude(digitID) ; // Calibration for an energy in GeV @@ -266,7 +296,45 @@ void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/) phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD); - Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); + // Hadronic Correction + double correction=0; + if (fApplyFractionHadronicCorrection) { + + TArrayS* matched = new TArrayS(); + GetTracksPointingToCell(matched, etaD, phiD,0.02); + + double ptot = 0; + for(int itrack = 0; itrack < matched->GetSize(); itrack++) + { + const short index = matched->At(itrack); + if (index>-1) + { + AliAODTrack *track = fAOD->GetTrack(index); + ptot += track->P(); + //printf("aod-track p=%f \n",track->P() ); + } + } + + correction = ptot * fFractionHadronicCorrection; + + if (ptot>0 && fDebug>1 ) { + printf("SUM of track momentum for this cell: p=%f \n", ptot); + printf("fractional correction=%f \n", fFractionHadronicCorrection); + printf("TOTAL correction=%f \n", correction ); + } + + delete matched; + + }//end hadronic correction + + double enerCorr = digitEn; + if (correction >= enerCorr) enerCorr = 0; + else enerCorr -= correction; + if (correction>0 && fDebug>1) { + printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr); + } + + Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID); if(uArray->GetUnitEnergy() == 0.) goodDigit++; @@ -317,6 +385,43 @@ void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/) } } +//____________________________________________________________________________ +void AliJetAODFillUnitArrayEMCalDigits::GetTracksPointingToCell(TArrayS* array,Double_t eta, Double_t phi, Double_t cut) +{ + + int size=0; + + for (Int_t itrk = 0; itrk < fAOD->GetNumberOfTracks() ; itrk++) { //track loop + + AliAODTrack * track = (AliAODTrack*) fAOD->GetTrack(itrk) ; + AliAODPid* pid = (AliAODPid*) track->GetDetPid(); + + if(pid) { + Double_t emcpos[3]; + pid->GetEMCALPosition(emcpos); + TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]); + + Double_t deta = tpos.Eta() - eta; + Double_t dphi = tpos.Phi() - phi; + + if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi(); + if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi(); + + Double_t res = sqrt(dphi*dphi + deta*deta); + + if (res< cut) { + //add this track-index + size++; + array->Set( size ); + array->AddAt( itrk, (size-1) ); + if(fDebug>1) printf("MTH:: track %d matched cell at eta=%f , phi=%f \n", itrk, eta, phi); + + } + } + } + +} + /* //____________________________________________________________________________ void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters() diff --git a/JETAN/AliJetAODFillUnitArrayEMCalDigits.h b/JETAN/AliJetAODFillUnitArrayEMCalDigits.h index eff0287165e..208cd4aac40 100755 --- a/JETAN/AliJetAODFillUnitArrayEMCalDigits.h +++ b/JETAN/AliJetAODFillUnitArrayEMCalDigits.h @@ -53,7 +53,9 @@ class AliJetAODFillUnitArrayEMCalDigits : public AliJetFillUnitArray Bool_t fApplyFractionHadronicCorrection; // Fraction hadronic correction flag Bool_t fFractionHadronicCorrection; // Fraction hadronic correction - + //Track-matching (mth) + void GetTracksPointingToCell(TArrayS *arr, Double_t eta, Double_t phi, Double_t res); + // geometry info AliESDCaloCluster *fClus; //! Int_t fNDigitEmcal; //! diff --git a/JETAN/AliJetAODFillUnitArrayTracks.cxx b/JETAN/AliJetAODFillUnitArrayTracks.cxx index c63816815ee..c4bd4edf71c 100644 --- a/JETAN/AliJetAODFillUnitArrayTracks.cxx +++ b/JETAN/AliJetAODFillUnitArrayTracks.cxx @@ -28,6 +28,7 @@ #include #include "AliJetUnitArray.h" +#include "AliJetHadronCorrectionv1.h" #include "AliJetAODFillUnitArrayTracks.h" // --- ROOT system --- @@ -109,7 +110,7 @@ AliJetAODFillUnitArrayTracks& AliJetAODFillUnitArrayTracks::operator=(const AliJ //____________________________________________________________________________ void AliJetAODFillUnitArrayTracks::InitParameters() { - fHadCorr = 0; // For hadron correction + // fHadCorr = 0; // For hadron correction fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, @@ -146,6 +147,7 @@ void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/) // Set parameters InitParameters(); + fRef->Clear(); // get number of tracks in event (for the loop) Int_t goodTrack = 0; @@ -157,8 +159,8 @@ void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/) TVector3 p3; nt = fAOD->GetNumberOfTracks(); - if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl; - + if(fDebug>1) cout << "Number of Tracks in AOD : " << nt << endl; + // temporary storage of signal and pt cut flag Int_t* sflag = new Int_t[nt]; Int_t* cflag = new Int_t[nt]; @@ -180,14 +182,8 @@ void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/) Int_t nTracksTpcCut = 0; Int_t nTracksTpcOnlyCut = 0; -//(not used ?) Int_t nElements = fTPCGrid->GetNEntries(); -//(not used ?) Int_t nElements2 = fEMCalGrid->GetNEntries(); fGrid = fTPCGrid->GetGridType(); -//(not used ?) Int_t nEntries = fUnitArray->GetEntries(); - -//(not used ?) Int_t nRefEnt = fRefArray->GetEntries(); - //loop over tracks nmax = nt; for (Int_t it = 0; it < nmax; it++) { @@ -201,8 +197,6 @@ void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/) p3.SetXYZ(mom[0],mom[1],mom[2]); pt = p3.Pt(); -//(not used ?) Float_t mass = 0.; - eta = p3.Eta(); phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi()); if (status == 0) continue; @@ -210,10 +204,13 @@ void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/) if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut + if (goodTrack == 0) new(fRef) TRefArray(TProcessID::GetProcessWithUID(track)); + sflag[goodTrack]=0; if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; cflag[goodTrack]=0; if (pt > ptMin) cflag[goodTrack]=1; // pt cut + fRef->Add(track); if(fGrid==0) { @@ -291,9 +288,7 @@ void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/) Int_t n1 = fGrid1->GetNEntries(); Int_t n2 = fGrid2->GetNEntries(); Int_t n3 = fGrid3->GetNEntries(); -//(not used ?) Int_t n4 = fGrid4->GetNEntries(); - // cout << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl; if(phi >= phimin0 && phi <= phimax0){ Int_t id0 = fGrid0->GetIndex(phi,eta)-1; AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements); @@ -518,15 +513,29 @@ void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/) } // Do Hadron Correction - // This is under construction !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - // Parametrization to be added - if (fApplyMIPCorrection != 0) + // For the moment I apply MIP correction if p >= 0.5 GeV/c + // and only for tracks inside EMCal acceptance + // End of if(fGrid==1) -> hadron correction for all tracks + if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5) { -// Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0); -// unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta))); + ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETE",1.); + + Double_t etaOut = 0.; + Double_t phiOut = 0.; + ((AliJetHadronCorrectionv1*)fHadCorr)->TrackPositionEMCal(track,etaOut,phiOut); + + // One has to make sure that the track hits the calorimeter + // We can then correct on average for these particles + if((etaOut >= fEtaMin && etaOut <= fEtaMax) && + (phiOut >= fPhiMin && phiOut <= fPhiMax))// && + { + Double_t hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0); + unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta))); + } } //end Hadron Correction loop - uArray->SetUnitEnergy(unitEnergy + pt); + if((unitEnergy + pt) > 0.) uArray->SetUnitEnergy(unitEnergy + pt); + else uArray->SetUnitEnergy(0.); // Put a pt cut flag if(uArray->GetUnitEnergy()Add(uArray); } } + + /* + // Do hadron correction + // In principle, the best thing to do as particles which + // eta,phi (at vertex) correspond to EMCal dead zone can deposit energy + // in the calorimeter as well as particles with eta,phi (at vertex) + // outside the EMCal acceptance + cout << "Hadronic correction for all tracks goes here" << endl; + + // For the moment I apply MIP correction if p >= 0.5 GeV/c + if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5) + { + ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETE",1.); + + Double_t etaOut = 0.; + Double_t phiOut = 0.; + Int_t towerID2 = 0; + Float_t unitEnergy = 0.; + ((AliJetHadronCorrectionv1*)fHadCorr)->TrackPositionEMCal(track,etaOut,phiOut); + + // One has to make sure that the track hits the calorimeter + // We can then correct on average for these particles + if((etaOut >= fEtaMin && etaOut <= fEtaMax) && + (phiOut >= fPhiMin && phiOut <= fPhiMax))// && + { + Int_t towerID = 0; + Int_t towerID2 = 0; + fGeom->GetAbsCellIdFromEtaPhi(etaOut,phiOut,towerID2); + if(towerID2==-1) continue; + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID2); + unitEnergy = uArray->GetUnitEnergy(); + Bool_t ok = kFALSE; + Bool_t ref = kFALSE; + if(unitEnergy==0.){ + nTracksEmcal++; + ok=kTRUE; + ref=kTRUE; + } + + Double_t hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0); + unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta))); + + // The energy in this unitarray can be negative + // If the Digits are then filled, make sure the energy sum + // is positive. If not, put the value to zero in AliJetAODUnitArrayEMCalDigit. + // If they are not filled, put this value to zero + uArray->SetUnitEnergy(unitEnergy); + + if(uArray->GetUnitEnergy()!=0. && ref){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + } + } + else cout << "Track out of EMCal acceptance" << endl; + + } //end Hadron Correction loop + */ + } // end fGrid==1 goodTrack++; @@ -609,6 +680,9 @@ void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/) cout << "goodTracks: " << goodTrack << endl; } + fSignalFlag.Set(goodTrack,sflag); + fCutFlag.Set(goodTrack,cflag); + delete[] sflag; delete[] cflag; diff --git a/JETAN/AliJetAODFillUnitArrayTracks.h b/JETAN/AliJetAODFillUnitArrayTracks.h index 67dd8fa0eb2..80609773149 100644 --- a/JETAN/AliJetAODFillUnitArrayTracks.h +++ b/JETAN/AliJetAODFillUnitArrayTracks.h @@ -46,7 +46,7 @@ class AliJetAODFillUnitArrayTracks : public AliJetFillUnitArray AliJetHadronCorrection *fHadCorr; // Pointer to Hadron Correction Object Bool_t fApplyMIPCorrection; // Apply MIP or not ? Exclusive with fApplyFractionHadronicCorrection - AliAODEvent *fAOD; // ESD + AliAODEvent *fAOD; // AOD output Event AliJetGrid *fGrid0; // Grid used for dead zones definition AliJetGrid *fGrid1; // Grid used for dead zones definition AliJetGrid *fGrid2; // Grid used for dead zones definition diff --git a/JETAN/AliJetAODReader.cxx b/JETAN/AliJetAODReader.cxx index fec35698e2e..3c10f28e57f 100644 --- a/JETAN/AliJetAODReader.cxx +++ b/JETAN/AliJetAODReader.cxx @@ -40,7 +40,7 @@ #include "AliJetDummyGeo.h" #include "AliJetAODFillUnitArrayTracks.h" #include "AliJetAODFillUnitArrayEMCalDigits.h" -#include "AliJetHadronCorrection.h" +#include "AliJetHadronCorrectionv1.h" #include "AliJetUnitArray.h" ClassImp(AliJetAODReader) @@ -424,7 +424,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); @@ -450,6 +450,7 @@ Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray) // clear momentum array ClearArray(); + fRef->Clear(); fDebug = fReaderHeader->GetDebug(); fOpt = fReaderHeader->GetDetector(); @@ -464,6 +465,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"); @@ -471,6 +475,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 @@ -528,13 +534,6 @@ 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 if(fDebug>1) printf("\n EMCal parameters initiated ! \n"); } diff --git a/JETAN/AliJetESDFillUnitArrayTracks.cxx b/JETAN/AliJetESDFillUnitArrayTracks.cxx index 050de40d869..05d41409b71 100644 --- a/JETAN/AliJetESDFillUnitArrayTracks.cxx +++ b/JETAN/AliJetESDFillUnitArrayTracks.cxx @@ -29,6 +29,7 @@ // --- AliRoot header files --- #include "AliJetUnitArray.h" +#include "AliJetHadronCorrectionv1.h" #include "AliJetESDFillUnitArrayTracks.h" // --- ROOT system --- @@ -112,7 +113,7 @@ AliJetESDFillUnitArrayTracks& AliJetESDFillUnitArrayTracks::operator=(const AliJ //____________________________________________________________________________ void AliJetESDFillUnitArrayTracks::InitParameters() { - fHadCorr = 0; // For hadron correction + // fHadCorr = 0; // For hadron correction fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL cout << "In AliJetESDFillUnitArrayTracks:InitParameters(), Ncells : " << fNumUnits << endl; @@ -149,6 +150,7 @@ void AliJetESDFillUnitArrayTracks::Exec(Option_t* const /*option*/) // Set parameters InitParameters(); + // fRef->Clear(); // get number of tracks in event (for the loop) Int_t goodTrack = 0; @@ -223,6 +225,7 @@ void AliJetESDFillUnitArrayTracks::Exec(Option_t* const /*option*/) if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; cflag[goodTrack]=0; if (pt > ptMin) cflag[goodTrack]=1; // pt cut + // fRef->Add(track); if(fGrid==0) { @@ -527,14 +530,30 @@ void AliJetESDFillUnitArrayTracks::Exec(Option_t* const /*option*/) // Do Hadron Correction // This is under construction !!!!!!!!!!!!!!!!!!!!!!! - // Parametrization to be added - if (fApplyMIPCorrection != 0) + // For the moment I apply MIP correction if p >= 0.5 GeV/c + if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5) { -// Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0); -// unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta))); + ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETE",1.); + + // Get track position at the outer part of the reconstruction ~ TRD + Double_t phiOut = track->GetOuterParam()->Phi(); + Double_t etaOut = track->GetOuterParam()->Eta(); + + // If the track in the outer part of the TPC/TDR ? is inside + // the calorimeter, it can deposit part of its energy + // We can then correct on average for these particles + if((etaOut >= fEtaMin && etaOut <= fEtaMax) && + (phiOut >= fPhiMin && phiOut <= fPhiMax))// && + { + Double_t hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0); + unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta))); + } } //end Hadron Correction loop - uArray->SetUnitEnergy(unitEnergy + pt); + cout << "unitEnergy + pt = " << unitEnergy << " + " << pt << " = " << unitEnergy + pt << endl; + + if((unitEnergy + pt) > 0.) uArray->SetUnitEnergy(unitEnergy + pt); + else uArray->SetUnitEnergy(0.); // Put a pt cut flag if(uArray->GetUnitEnergy()SetEMCalGrid(fEmcalGrid); fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection); fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection); - fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection); + fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection); // Add the task to global task fFillUnitArray->Add(fFillUAFromTracks); fFillUnitArray->Add(fFillUAFromEMCalDigits); @@ -430,7 +431,7 @@ void AliJetESDReader::InitParameters() { // Initialise parameters fOpt = fReaderHeader->GetDetector(); - fHadCorr = 0; // For hadron correction + // fHadCorr = 0; // For hadron correction if(fEFlag==kFALSE){ if(fOpt==0 || fOpt==1) fECorrection = 0; // For electron correction diff --git a/JETAN/AliJetFillUnitArray.cxx b/JETAN/AliJetFillUnitArray.cxx index e2a91abf9a4..628808065a4 100644 --- a/JETAN/AliJetFillUnitArray.cxx +++ b/JETAN/AliJetFillUnitArray.cxx @@ -53,6 +53,9 @@ AliJetFillUnitArray::AliJetFillUnitArray() fMomentumArray(0x0), fUnitArray(0x0), fRefArray(0x0), + fRef(0x0), + fSignalFlag(0), + fCutFlag(0), fProcId(kFALSE), fTPCGrid(0x0), fEMCalGrid(0x0), diff --git a/JETAN/AliJetFillUnitArray.h b/JETAN/AliJetFillUnitArray.h index 1c190b57ed3..8dabec5f23c 100644 --- a/JETAN/AliJetFillUnitArray.h +++ b/JETAN/AliJetFillUnitArray.h @@ -15,6 +15,7 @@ #include #include +#include #include #include @@ -46,6 +47,9 @@ class AliJetFillUnitArray : public TTask virtual void SetMomentumArray(TClonesArray* const momentumArray) {fMomentumArray = momentumArray;} virtual void SetUnitArray(TClonesArray* const unitArray) {fUnitArray = unitArray;} virtual void SetRefArray(TRefArray* const refArray) {fRefArray = refArray;} + virtual void SetReferences(TRefArray *refArray) {fRef = refArray;} + virtual void SetSignalFlag(TArrayI sflag) {fSignalFlag = sflag;} + virtual void SetCutFlag(TArrayI cflag) {fCutFlag = cflag;} virtual void SetTPCGrid(AliJetGrid* const grid) {fTPCGrid = grid;} virtual void SetEMCalGrid(AliJetGrid* const grid) {fEMCalGrid = grid;} virtual void SetProcId(Bool_t id) {fProcId = id;} @@ -73,6 +77,8 @@ class AliJetFillUnitArray : public TTask virtual Int_t GetMult() const {return 0;} // To be checked virtual Int_t GetMultCut() const {return 0;} // To be checked virtual Bool_t GetProcId() const {return fProcId;} + virtual TArrayI GetSignalFlag() const {return fSignalFlag;} + virtual TArrayI GetCutFlag() const {return fCutFlag;} // Other virtual void Exec(Option_t* const /*option*/) {;} @@ -90,6 +96,11 @@ class AliJetFillUnitArray : public TTask TClonesArray *fMomentumArray; // MomentumArray TClonesArray *fUnitArray; // UnitArray TRefArray *fRefArray; // UnitArray + TRefArray *fRef; // ref Array to aod tracks + TArrayI fSignalFlag; // to flag if a particle comes from pythia or + // from the underlying event + TArrayI fCutFlag; // to flag if a particle passed the pt cut or not + Bool_t fProcId; // Bool_t for TProcessID synchronization AliJetGrid *fTPCGrid; // Define filled grid AliJetGrid *fEMCalGrid; // Define filled grid diff --git a/JETAN/AliJetHadronCorrectionv0.cxx b/JETAN/AliJetHadronCorrectionv0.cxx index 6fe34e88e99..5d65d6b98ba 100644 --- a/JETAN/AliJetHadronCorrectionv0.cxx +++ b/JETAN/AliJetHadronCorrectionv0.cxx @@ -26,10 +26,10 @@ #include "AliJetHadronCorrectionv0.h" const Int_t maxVariant = 8; // size eta grid -const Int_t nVec = 10; // size momentum grid -const Int_t nPol = 4; // number coefficients of polinom +const Int_t nVec = 11; // size momentum grid +const Int_t nPol = 4; // number coefficients of polynom static Double_t etaGrid[maxVariant]={ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.67}; -static Double_t pGrid[nVec]={ 0.2, 0.5, 1.0, 2.0, 3.0, 5.0,10.0,15.0,25.0,40.0}; +static Double_t pGrid[nVec]={ 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 30.0, 50.0, 70.0, 100.0}; // c[][] - first index for eta, second for momentum static Double_t c[maxVariant][nPol] ={ {1.305705e-01, 3.725653e-01, -1.219962e-02, 1.806235e-04}, @@ -59,7 +59,7 @@ AliJetHadronCorrectionv0::Instance() } Double_t -AliJetHadronCorrectionv0::GetEnergy(Double_t pmom, Double_t eta, Int_t /*gid*/) +AliJetHadronCorrectionv0::GetEnergy(Double_t pmom, Double_t eta, Int_t gid) { Int_t iEta=0; // index Double_t etaw = TMath::Abs(eta); diff --git a/JETAN/AliJetHadronCorrectionv0.h b/JETAN/AliJetHadronCorrectionv0.h index 998fbc16049..d8a3001d1bd 100644 --- a/JETAN/AliJetHadronCorrectionv0.h +++ b/JETAN/AliJetHadronCorrectionv0.h @@ -22,8 +22,7 @@ class AliJetHadronCorrectionv0: public AliJetHadronCorrection { public: static AliJetHadronCorrectionv0* Instance(); virtual Double_t GetEnergy(Double_t pmom, Double_t eta, Int_t gid); - Double_t GetEnergy(Double_t pmom, Double_t eta) - {return GetEnergy(pmom,eta,7);} + Double_t GetEnergy(Double_t pmom, Double_t eta){return GetEnergy(pmom,eta,7);} virtual ~AliJetHadronCorrectionv0() {} diff --git a/JETAN/AliJetHadronCorrectionv1.cxx b/JETAN/AliJetHadronCorrectionv1.cxx index 920c8b95496..b9341f5b91f 100644 --- a/JETAN/AliJetHadronCorrectionv1.cxx +++ b/JETAN/AliJetHadronCorrectionv1.cxx @@ -15,7 +15,7 @@ //=============================================================== // To be modified for hadron correction using particle ID and track merging -// Author : magali.estienne@ires.in2p3.fr +// Author : magali.estienne@subatech.in2p3.fr //=============================================================== // Author : Mark Horner (LBL/UCT) @@ -24,37 +24,55 @@ #include "TMath.h" // --- AliRoot header files --- -//#include "../EMCAL/AliJetGeometry.h" +#include "AliAODTrack.h" #include "AliJetDummyGeo.h" #include "AliJetHadronCorrectionv1.h" +static Double_t etaGrid[HCPARAMETERS]={ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.67}; + ClassImp(AliJetHadronCorrectionv1) Double_t AliJetHadronCorrectionv1::fgParLookup[HCPARAMETERS][HCPARAMETERSETS] = { - {-2.82271e-4 , -2.39954e-4}, - {2.50796e-2 , 2.07172e-2}, - {1.02861e-3 , 1.48576e-3}, - {2.11539e-2 , -1.38473-2}, - {2.27003e-2 , 2.78252e-2}, - {1.65078e-6 , 1.51821e-6} + + {7.52624e-2 , 9.80119e-2 , 1.29086e-2}, + {-2.07449e-2 , 2.22951e-2 , 8.31971e-2}, + {-3.64861e-4 , -3.65680e-03 , -3.07148e-02}, + {5.35252e-3 , 5.12234e-03 , 2.25559e-02}, + {1.27106e-1 , 1.25248e-01 , 1.06352e-01}, + {-6.72909e-2 , 9.78677e-01 , -1.35909e-01}, + {-3.80660e-5 , -2.77522e-05 , -2.05218e-05}, + {1.54855e-7 , 1.00604e-07 , 1.02481e-07} }; AliJetHadronCorrectionv1* AliJetHadronCorrectionv1::fgHadrCorr = 0; -void AliJetHadronCorrectionv1::SetGeometry(AliJetDummyGeo *geometry) +AliJetHadronCorrectionv1::AliJetHadronCorrectionv1(const char *name,const char *title) + :AliJetHadronCorrection(name, title), + fSamplingFraction(0) +{ + fgHadrCorr = this; +} + +AliJetHadronCorrectionv1* +AliJetHadronCorrectionv1::Instance() +{ + // return pointer to global instance. Instantiate if needed + if (! fgHadrCorr) fgHadrCorr = new AliJetHadronCorrectionv1(); + return fgHadrCorr; +} + +void AliJetHadronCorrectionv1::SetGeometry2(AliJetDummyGeo *geometry) { // Initialise EMCAL geometry if (!geometry) { SetParameters(); fSamplingFraction = geometry->GetSampling(); - cout<<"Setting the sampling fraction to :"<GetName()); fSamplingFraction = geometry->GetSampling(); - cout<<"Setting the sampling fraction to :"<= etaGrid[1] && etai <= etaGrid[HCPARAMETERS-2]) { + for (Int_t i=0;i<8;i++) + { + fPar[i] = fgParLookup[i][0]; + // cout << "fPar[" << i << "]: " << fPar[i] << endl; + } + } else { + for (Int_t i=0;i<8;i++) + { + fPar[i] = fgParLookup[i][2]; + // cout << "fPar[" << i << "]: " << fPar[i] << endl; + } + + } + + Double_t value = fPar[5]*pow(etai,3) + + fPar[0]*pow(etai,2) + + fPar[1]*etai + + fPar[2]*etai*pmom + + fPar[3]*pmom + + fPar[4] + + fPar[6]*pow(pmom,2) + + fPar[7]*pow(pmom,3); + return fSamplingFraction*value; } + +void AliJetHadronCorrectionv1::TrackPositionEMCal(AliAODTrack* track,Double_t &eta, Double_t &phi) +{ + AliAODPid* pid = (AliAODPid*) track->GetDetPid(); + + if(pid) { + Double_t emcpos[3]; + pid->GetEMCALPosition(emcpos); + TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]); + + eta = tpos.Eta(); + phi = tpos.Phi(); + + } + +} diff --git a/JETAN/AliJetHadronCorrectionv1.h b/JETAN/AliJetHadronCorrectionv1.h index 448c42bd1fb..b837319e20e 100644 --- a/JETAN/AliJetHadronCorrectionv1.h +++ b/JETAN/AliJetHadronCorrectionv1.h @@ -8,36 +8,37 @@ // #include "AliJetHadronCorrection.h" -#define HCPARAMETERS 6 -#define HCPARAMETERSETS 2 +#define HCPARAMETERS 8 +#define HCPARAMETERSETS 3 class AliJetDummyGeo; - class AliJetHadronCorrectionv1 : public AliJetHadronCorrection { public: - AliJetHadronCorrectionv1(); - static AliJetHadronCorrectionv1* Instance(); - virtual Double_t GetEnergy(Double_t pmom, Double_t eta, Int_t gid); - Double_t GetEnergy(Double_t pmom, Double_t eta) {return GetEnergy(pmom,eta,7);} - - void SetGeometry(TString name, Double_t fs = 1.); - virtual ~AliJetHadronCorrectionv1() {} - + static AliJetHadronCorrectionv1* Instance(); + virtual ~AliJetHadronCorrectionv1() {} + + virtual Double_t GetEnergy(Double_t pmom, Double_t eta, Int_t gid); + Double_t GetEnergy(Double_t pmom, Double_t eta){return GetEnergy(pmom,eta,7);} + + void SetGeometry(TString name, Double_t fs = 1.); + void SetGeometry2(AliJetDummyGeo *geometry); + void TrackPositionEMCal(AliAODTrack* track,Double_t &eta, Double_t &phi); + protected: - AliJetHadronCorrectionv1(const char *name, const char *title); - void SetGeometry(AliJetDummyGeo *geometry); - + AliJetHadronCorrectionv1() {;} + AliJetHadronCorrectionv1(const char *name, const char *title); + private: - void SetParameters(TString name = "") {Warning("SetParameter","Dummy method with argument %s",name.Data());} - - static AliJetHadronCorrectionv1* fgHadrCorr; // Pointer to global instance (singleton) - static Double_t fgParLookup[HCPARAMETERS][HCPARAMETERSETS]; // Global array with parameters for hadronic response - Double_t fPar[6]; - Float_t fSamplingFraction; // Sampling fraction + void SetParameters(TString name = "") {Warning("SetParameter","Dummy method with argument %s",name.Data());} + + static AliJetHadronCorrectionv1* fgHadrCorr; // Pointer to global instance (singleton) + static Double_t fgParLookup[HCPARAMETERS][HCPARAMETERSETS]; // Global array with parameters for hadronic response + Double_t fPar[8]; + Float_t fSamplingFraction; // Sampling fraction - ClassDef(AliJetHadronCorrectionv1,2) // Hadron correction for EMC (version for MDC) + ClassDef(AliJetHadronCorrectionv1,2) // Hadron correction for EMC (version for MDC) }; #endif // ALIJETHADRONCORRECTIONV1_H diff --git a/JETAN/AliUA1JetFinderV2.cxx b/JETAN/AliUA1JetFinderV2.cxx index c401860be26..074dd4b4366 100644 --- a/JETAN/AliUA1JetFinderV2.cxx +++ b/JETAN/AliUA1JetFinderV2.cxx @@ -613,6 +613,14 @@ void AliUA1JetFinderV2::FindJets() multJetOk[p] = multJet[idx[p]]; } + TRefArray *refs = 0; + Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader"); + if (fromAod) refs = fReader->GetReferences(); + Int_t nTracks = 0; + if (fromAod) nTracks = ((TRefArray*)refs)->GetEntries(); + Int_t* trackinjet = new Int_t[nTracks]; + for(Int_t it=0; it (header->GetJetEtaMax())) || @@ -627,6 +635,22 @@ void AliUA1JetFinderV2::FindJets() AliAODJet jet(px, py, pz, en); jet.Print(""); + if (fromAod){ + for(Int_t jpart = 0; jpart < nTracks; jpart++) { // loop for all particles in array + Float_t deta = ((AliAODTrack*)refs->At(jpart))->Eta() - etaJetOk[kj]; + Float_t dphi = ((AliAODTrack*)refs->At(jpart))->Phi() - phiJetOk[kj]; + if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); + if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= header->GetRadius() && fReader->GetCutFlag(jpart) == 1) // particles inside this cone + if(trackinjet[jpart]==-1) trackinjet[jpart] = kj; + else if(fDebug>10) printf("The track already belongs to jet %d \n",trackinjet[jpart]); + if(trackinjet[jpart]==kj) + jet.AddTrack(refs->At(jpart)); // check if the particle belongs to the jet and add the ref + } + } + AddJet(jet); idxjets[nselectj] = kj; @@ -701,6 +725,7 @@ void AliUA1JetFinderV2::FindJets() delete ncellsJetOk; delete multJetOk; //-------------------------- + delete trackinjet; delete idxjets; delete percentage; delete ncells; -- 2.43.0