fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
fPosAlgo(kUnchanged),fW0(4.),
fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
- fRemoveBadChannels(kFALSE), fEMCALBadChannelMap(),
+ fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
fMatchedClusterIndex(0x0), fResidualZ(0x0), fResidualR(0x0), fCutR(20), fCutZ(20),
fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
: TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
- fRemoveBadChannels(reco.fRemoveBadChannels),fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
+ fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
+ fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
fRecalibration = reco.fRecalibration;
fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
fRemoveBadChannels = reco.fRemoveBadChannels;
+ fRecalDistToBadChannels= reco.fRecalDistToBadChannels;
fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
//If the distance to the border is 0 or negative just exit accept all clusters
if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
- Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
- GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi);
+ Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
+ Bool_t shared = kFALSE;
+ GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
switch (fNonLinearityFunction) {
case kPi0MC:
+ {
//Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
- //Double_t par0 = 1.001;
- //Double_t par1 = -0.01264;
- //Double_t par2 = -0.03632;
- //Double_t par3 = 0.1798;
- //Double_t par4 = -0.522;
+ //Double_t fNonLinearityParams[0] = 1.001;
+ //Double_t fNonLinearityParams[1] = -0.01264;
+ //Double_t fNonLinearityParams[2] = -0.03632;
+ //Double_t fNonLinearityParams[3] = 0.1798;
+ //Double_t fNonLinearityParams[4] = -0.522;
energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
break;
+ }
case kPi0GammaGamma:
-
+ {
//Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
- //Double_t par0 = 0.1457;
- //Double_t par1 = -0.02024;
- //Double_t par2 = 1.046;
+ //Double_t fNonLinearityParams[0] = 0.1457;
+ //Double_t fNonLinearityParams[1] = -0.02024;
+ //Double_t fNonLinearityParams[2] = 1.046;
energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
break;
+ }
case kPi0GammaConversion:
-
+ {
//Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
- //Double_t C = 0.139393/0.1349766;
- //Double_t a = 0.0566186;
- //Double_t b = 0.982133;
+ //fNonLinearityParams[0] = 0.139393/0.1349766;
+ //fNonLinearityParams[1] = 0.0566186;
+ //fNonLinearityParams[2] = 0.982133;
energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
break;
+ }
+
+ case kBeamTest:
+ {
+ //From beam test, Alexei's results, for different ZS thresholds
+ // th=30 MeV; th = 45 MeV; th = 75 MeV
+ //fNonLinearityParams[0] = 0.107; 1.003; 1.002
+ //fNonLinearityParams[1] = 0.894; 0.719; 0.797
+ //fNonLinearityParams[2] = 0.246; 0.334; 0.358
+ energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
+
+ break;
+ }
case kNoCorrection:
AliDebug(2,"No correction on the energy\n");
//Calculate shower depth for a given cluster energy and particle type
// parameters
- Float_t x0 = 1.23;
+ Float_t x0 = 1.31;
Float_t ecr = 8;
Float_t depth = 0;
}
//__________________________________________________
-void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
+void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
+ Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
{
//For a given CaloCluster gets the absId of the cell
//with maximum energy deposit.
Int_t iTower = -1;
Int_t iIphi = -1;
Int_t iIeta = -1;
+ Int_t iSupMod0= -1;
//printf("---Max?\n");
for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
cellAbsId = clu->GetCellAbsId(iDig);
fraction = clu->GetCellAmplitudeFraction(iDig);
//printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
+ geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
+ geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
+ if(iDig==0) iSupMod0=iSupMod;
+ else if(iSupMod0!=iSupMod) {
+ shared = kTRUE;
+ //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
+ }
if(IsRecalibrationOn()) {
- geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
- geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
}
eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
- fEMCALRecalibrationFactors = new TObjArray(12);
+ fEMCALRecalibrationFactors = new TObjArray(10);
for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
//Init the histograms with 1
for (Int_t sm = 0; sm < 12; sm++) {
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
- fEMCALBadChannelMap = new TObjArray(12);
+ fEMCALBadChannelMap = new TObjArray(10);
//TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
- for (int i = 0; i < 12; i++) {
+ for (int i = 0; i < 10; i++) {
fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
}
Float_t weight = 0., totalWeight=0.;
Float_t newPos[3] = {0,0,0};
Double_t pLocal[3], pGlobal[3];
-
+ Bool_t shared = kFALSE;
+
Float_t clEnergy = clu->E(); //Energy already recalibrated previously
- GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi);
+ GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
//printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
Int_t iIphi = -1, iIeta = -1;
Int_t iSupMod = -1, iSupModMax = -1;
Int_t iphi = -1, ieta =-1;
-
+ Bool_t shared = kFALSE;
+
Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
- GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi);
+ GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
}
+//____________________________________________________________________________
+void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
+
+ //re-evaluate distance to bad channel with updated bad map
+
+ if(!fRecalDistToBadChannels) return;
+
+ //Get channels map of the supermodule where the cluster is.
+ Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
+ Bool_t shared = kFALSE;
+ GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
+ TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
+
+ Int_t dRrow, dRcol;
+ Float_t minDist = 10000.;
+ Float_t dist = 0.;
+
+ //Loop on tower status map
+ for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
+ for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
+ //Check if tower is bad.
+ if(hMap->GetBinContent(icol,irow)==0) continue;
+ //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
+ // iSupMod,icol, irow, icolM,irowM);
+
+ dRrow=TMath::Abs(irowM-irow);
+ dRcol=TMath::Abs(icolM-icol);
+ dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
+ if(dist < minDist){
+ //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
+ minDist = dist;
+ }
+
+ }
+ }
+
+ //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
+ if (shared) {
+ TH2D* hMap2 = 0;
+ Int_t iSupMod2 = -1;
+
+ //The only possible combinations are (0,1), (2,3) ... (8,9)
+ if(iSupMod%2) iSupMod2 = iSupMod-1;
+ else iSupMod2 = iSupMod+1;
+ hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
+
+ //Loop on tower status map of second super module
+ for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
+ for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
+ //Check if tower is bad.
+ if(hMap2->GetBinContent(icol,irow)==0) continue;
+ //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels(shared) - \n \t Bad channel in SM %d, col %d, row %d \n \t Cluster max in SM %d, col %d, row %d\n",
+ // iSupMod2,icol, irow,iSupMod,icolM,irowM);
+
+ dRrow=TMath::Abs(irow-irowM);
+
+ if(iSupMod%2) {
+ dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
+ }
+ else {
+ dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
+ }
+
+ dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
+ if(dist < minDist) minDist = dist;
+
+ }
+ }
+
+ }// shared cluster in 2 SuperModules
+
+ AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
+ cluster->SetDistanceToBadChannel(minDist);
+
+}
+
//____________________________________________________________________________
void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
//Get the residuals dR and dZ for this cluster
//It only works with ESDs, not AODs
- if( FindMatchedPos(index)==-1 )
+ if( FindMatchedPos(index) >= 999 )
{
AliDebug(2,"No matched tracks found!\n");
dR=999.;
{
//Given a cluster index as in AliESDEvent::GetCaloCluster(index)
//Returns if cluster has a match
- if(FindMatchedPos(index)>-1) return kTRUE;
+ if(FindMatchedPos(index) < 999)
+ return kTRUE;
else
return kFALSE;
}
//__________________________________________________
-Int_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
+UInt_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
{
//Given a cluster index as in AliESDEvent::GetCaloCluster(index)
//Returns the position of the match in the fMatchedClusterIndex array
Float_t tmpR = fCutR;
- Int_t pos=-1;
+ UInt_t pos = 999;
for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
{
pos=i;
tmpR=fResidualR->At(i);
}
- //printf("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i));
+ AliDebug(3,Form("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
}
return pos;
}
chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
if (nClustersTPC!=0)
chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
-
+
+
+ //DCA cuts
+ Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
+ //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
+ SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
+
+
Float_t b[2];
Float_t bCov[3];
esdTrack->GetImpactParameters(b,bCov);
if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
cuts[9] = kTRUE;
+ //Require at least one SPD point + anything else in ITS
+ if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
+ cuts[10] = kTRUE;
+
Bool_t cut=kFALSE;
for (Int_t i=0; i<kNCuts; i++)
if (cuts[i]) {cut = kTRUE;}
void AliEMCALRecoUtils::InitTrackCuts()
{
//Intilize the track cut criteria
- //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
+ //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
//Also you can customize the cuts using the setters
-
- SetMinNClustersTPC(50);
+
+ //TPC
+ SetMinNClustersTPC(70);
SetMaxChi2PerClusterTPC(4);
SetAcceptKinkDaughters(kFALSE);
- SetMaxDCAToVertexZ(3.2);
- SetMaxDCAToVertexXY(2.4);
- SetDCAToVertex2D(kTRUE);
+ SetRequireTPCRefit(kTRUE);
+
+ //ITS
+ SetRequireITSRefit(kTRUE);
+ SetMaxDCAToVertexZ(2);
+ SetDCAToVertex2D(kFALSE);
}
//__________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff