AliPhysics  e6c8d43 (e6c8d43)
AliEMCALRecoUtils.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // ROOT includes
17 #include <TH2F.h>
18 #include <TArrayI.h>
19 #include <TArrayF.h>
20 #include <TObjArray.h>
21 
22 // STEER includes
23 #include "AliVCluster.h"
24 #include "AliVCaloCells.h"
25 #include "AliLog.h"
26 #include "AliPID.h"
27 #include "AliESDEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliESDtrack.h"
30 #include "AliAODTrack.h"
31 #include "AliExternalTrackParam.h"
32 #include "AliESDfriendTrack.h"
33 #include "AliTrackerBase.h"
34 #include "AliMCEvent.h"
35 
36 // EMCAL includes
37 #include "AliEMCALRecoUtils.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliTrackerBase.h"
40 #include "AliEMCALPIDUtils.h"
41 
43 ClassImp(AliEMCALRecoUtils) ;
45 
51 //_____________________________________
53  fParticleType(0), fPosAlgo(0),
54  fW0(0), fShowerShapeCellLocationType(0),
55  fNonLinearityFunction(0), fNonLinearThreshold(0),
56  fSmearClusterEnergy(kFALSE), fRandom(),
57  fCellsRecalibrated(kFALSE), fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
58  fConstantTimeShift(0), fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(), fLowGain(kFALSE),
59  fUseL1PhaseInTimeRecalibration(kFALSE), fEMCALL1PhaseInTimeRecalibration(),
60  fUseRunCorrectionFactors(kFALSE),
61  fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
62  fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
63  fRejectExoticCluster(kFALSE), fRejectExoticCells(kFALSE),
64  fExoticCellFraction(0), fExoticCellDiffTime(0), fExoticCellMinAmplitude(0),
65  fPIDUtils(), fAODFilterMask(0),
66  fAODHybridTracks(0), fAODTPCOnlyTracks(0),
67  fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
68  fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kFALSE), fCutEtaPhiSeparate(kFALSE),
69  fCutR(0), fCutEta(0), fCutPhi(0),
70  fClusterWindow(0), fMass(0),
71  fStepSurface(0), fStepCluster(0),
72  fITSTrackSA(kFALSE), fEMCalSurfaceDistance(440.),
73  fTrackCutsType(0), fCutMinTrackPt(0), fCutMinNClusterTPC(0),
74  fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
75  fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
76  fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0), fCutDCAToVertex2D(kFALSE),
77  fCutRequireITSStandAlone(kFALSE), fCutRequireITSpureSA(kFALSE),
78  fNMCGenerToAccept(0), fMCGenerToAcceptForTrack(1)
79 {
80  // Init parameters
82 
83  for(Int_t j = 0; j < 5; j++) fMCGenerToAccept[j] = "";
84 
85  // Track matching arrays init
88  fResidualPhi = new TArrayF();
89  fResidualEta = new TArrayF();
90  fPIDUtils = new AliEMCALPIDUtils();
91 }
92 
93 //
94 // Copy constructor.
95 //
96 //______________________________________________________________________
98 : AliEMCALRecoUtilsBase(reco),
107  fLowGain(reco.fLowGain),
112  fEMCALBadChannelMap(NULL),
121  fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0),
122  fResidualPhi( reco.fResidualPhi? new TArrayF(*reco.fResidualPhi):0x0),
124  fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
136 {
137  for (Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ;
138  fMisalTransShift[i] = reco.fMisalTransShift[i] ; }
139  for (Int_t i = 0; i < 10 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
140  for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
141  for (Int_t j = 0; j < 5 ; j++) { fMCGenerToAccept[j] = reco.fMCGenerToAccept[j] ; }
142 
143  if(reco.fEMCALBadChannelMap) {
144  // Copy constructor - not taking ownership over calibration histograms
145  fEMCALBadChannelMap = new TObjArray(reco.fEMCALBadChannelMap->GetEntries());
146  fEMCALBadChannelMap->SetOwner(false);
147  for(int ism = 0; ism < reco.fEMCALBadChannelMap->GetEntries(); ism++) fEMCALBadChannelMap->AddAt(reco.fEMCALBadChannelMap->At(ism), ism);
148  }
149 
150  if(reco.fEMCALRecalibrationFactors) {
151  // Copy constructor - not taking ownership over calibration histograms
153  fEMCALRecalibrationFactors->SetOwner(false);
154  for(int ism = 0; ism < reco.fEMCALRecalibrationFactors->GetEntries(); ism++) fEMCALRecalibrationFactors->AddAt(reco.fEMCALRecalibrationFactors->At(ism), ism);
155  }
156 
158  // Copy constructor - not taking ownership over calibration histograms
160  fEMCALTimeRecalibrationFactors->SetOwner(false);
161  for(int ism = 0; ism < reco.fEMCALTimeRecalibrationFactors->GetEntries(); ism++) fEMCALTimeRecalibrationFactors->AddAt(reco.fEMCALTimeRecalibrationFactors->At(ism), ism);
162  }
163 }
164 
168 //______________________________________________________________________
170 {
171  if (this == &reco)return *this;
172  ((TNamed *)this)->operator=(reco);
173 
174  for (Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ;
175  fMisalRotShift[i] = reco.fMisalRotShift[i] ; }
176  for (Int_t i = 0; i < 10 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
177  for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
178 
180  fPosAlgo = reco.fPosAlgo;
181  fW0 = reco.fW0;
183 
187 
190 
193  fLowGain = reco.fLowGain;
194 
197 
199 
202 
205 
211 
212  fPIDUtils = reco.fPIDUtils;
213 
217 
220  fCutR = reco.fCutR;
221  fCutEta = reco.fCutEta;
222  fCutPhi = reco.fCutPhi;
224  fMass = reco.fMass;
225  fStepSurface = reco.fStepSurface;
226  fStepCluster = reco.fStepCluster;
227  fITSTrackSA = reco.fITSTrackSA;
229 
244 
247  for (Int_t j = 0; j < 5 ; j++)
248  fMCGenerToAccept[j] = reco.fMCGenerToAccept[j];
249 
250  //
251  // Assign or copy construct the different TArrays
252  //
253  if (reco.fResidualEta)
254  {
255  if (fResidualEta)
256  *fResidualEta = *reco.fResidualEta;
257  else
258  fResidualEta = new TArrayF(*reco.fResidualEta);
259  }
260  else
261  {
262  if (fResidualEta) delete fResidualEta;
263  fResidualEta = 0;
264  }
265 
266  if (reco.fResidualPhi)
267  {
268  if (fResidualPhi)
269  *fResidualPhi = *reco.fResidualPhi;
270  else
271  fResidualPhi = new TArrayF(*reco.fResidualPhi);
272  }
273  else
274  {
275  if (fResidualPhi) delete fResidualPhi;
276  fResidualPhi = 0;
277  }
278 
279  if (reco.fMatchedTrackIndex)
280  {
281  if (fMatchedTrackIndex)
283  else
285  }
286  else
287  {
289  fMatchedTrackIndex = 0;
290  }
291 
292  if (reco.fMatchedClusterIndex)
293  {
296  else
298  }
299  else
300  {
303  }
304 
306  if(reco.fEMCALBadChannelMap) {
307  // Copy constructor - not taking ownership over calibration histograms
308  fEMCALBadChannelMap = new TObjArray(reco.fEMCALBadChannelMap->GetEntries());
309  fEMCALBadChannelMap->SetOwner(false);
310  for(int ism = 0; ism < reco.fEMCALBadChannelMap->GetEntries(); ism++) fEMCALBadChannelMap->AddAt(reco.fEMCALBadChannelMap->At(ism), ism);
311  }
312 
314  if(reco.fEMCALRecalibrationFactors) {
315  // Copy constructor - not taking ownership over calibration histograms
317  fEMCALRecalibrationFactors->SetOwner(false);
318  for(int ism = 0; ism < reco.fEMCALRecalibrationFactors->GetEntries(); ism++) fEMCALRecalibrationFactors->AddAt(reco.fEMCALRecalibrationFactors->At(ism), ism);
319  }
320 
323  // Copy constructor - not taking ownership over calibration histograms
325  fEMCALTimeRecalibrationFactors->SetOwner(false);
326  for(int ism = 0; ism < reco.fEMCALTimeRecalibrationFactors->GetEntries(); ism++) fEMCALTimeRecalibrationFactors->AddAt(reco.fEMCALTimeRecalibrationFactors->At(ism), ism);
327  }
328 
331  // Copy constructor - not taking ownership over calibration histograms
333  fEMCALL1PhaseInTimeRecalibration->SetOwner(false);
334  for(int ism = 0; ism < reco.fEMCALL1PhaseInTimeRecalibration->GetEntries(); ism++) fEMCALL1PhaseInTimeRecalibration->AddAt(reco.fEMCALL1PhaseInTimeRecalibration->At(ism), ism);
335  }
336 
337  return *this;
338 }
339 
343 //_____________________________________
345 {
347  {
349  }
350 
352  {
354  }
355 
357  {
359  }
360 
361  if (fEMCALBadChannelMap)
362  {
363  delete fEMCALBadChannelMap;
364  }
365 
366  delete fMatchedTrackIndex ;
367  delete fMatchedClusterIndex ;
368  delete fResidualEta ;
369  delete fResidualPhi ;
370  delete fPIDUtils ;
371 
372  InitTrackCuts();
373 }
374 
387 //_______________________________________________________________________________
389  Float_t & amp, Double_t & time,
390  AliVCaloCells* cells)
391 {
392  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
393 
394  if(!geom)
395  {
396  AliError("No instance of the geometry is available");
397  return kFALSE;
398  }
399 
400  if ( absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules() )
401  return kFALSE;
402 
403  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
404 
405  if (!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
406  {
407  // cell absID does not exist
408  amp=0; time = 1.e9;
409  return kFALSE;
410  }
411 
412  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
413 
414  // Do not include bad channels found in analysis,
415  if (IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi))
416  return kFALSE;
417 
418  //Recalibrate energy
419  amp = cells->GetCellAmplitude(absID);
421  amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
422 
423  // Recalibrate time
424  time = cells->GetCellTime(absID);
425  time-=fConstantTimeShift*1e-9; // only in case of old Run1 simulation
426  Bool_t isLowGain = !(cells->GetCellHighGain(absID));//HG = false -> LG = true
427 
428  RecalibrateCellTime(absID,bc,time,isLowGain);
429 
430  //Recalibrate time with L1 phase
431  RecalibrateCellTimeL1Phase(imod, bc, time);
432 
433  return kTRUE;
434 }
435 
445 //_____________________________________________________________________________
447  const AliVCluster* cluster,
448  AliVCaloCells* cells)
449 {
450  if (!cluster)
451  {
452  AliInfo("Cluster pointer null!");
453  return kFALSE;
454  }
455 
456  // If the distance to the border is 0 or negative just exit accept all clusters
457  if (cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 )
458  return kTRUE;
459 
460  Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
461  Bool_t shared = kFALSE;
462  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
463 
464  AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
465  absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
466 
467  if (absIdMax==-1) return kFALSE;
468 
469  // Check if the cell is close to the borders:
470  Bool_t okrow = kFALSE;
471  Bool_t okcol = kFALSE;
472 
473  if (iSM < 0 || iphi < 0 || ieta < 0 )
474  {
475  AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
476  iSM,ieta,iphi));
477  return kFALSE; // trick coverity
478  }
479 
480  // Check rows/phi
481  Int_t iPhiLast = 24;
482  if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_Half ) iPhiLast /= 2;
483  else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_3rd ) iPhiLast /= 3;// 1/3 sm case
484  else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Ext ) iPhiLast /= 3;// 1/3 sm case
485 
486  if(iphi >= fNCellsFromEMCALBorder && iphi < iPhiLast - fNCellsFromEMCALBorder) okrow = kTRUE;
487 
488  // Check columns/eta
489  Int_t iEtaLast = 48;
490  if ( !fNoEMCALBorderAtEta0 || geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Standard )
491  {
492  // consider inner border
493  if ( geom->IsDCALSM(iSM) ) iEtaLast = iEtaLast*2/3;
494 
495  if ( ieta > fNCellsFromEMCALBorder && ieta < iEtaLast-fNCellsFromEMCALBorder ) okcol = kTRUE;
496  }
497  else
498  {
499  if (iSM%2==0)
500  {
501  if (ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
502  }
503  else
504  {
505  if (ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE;
506  }
507  }//eta 0 not checked
508 
509  AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
510  fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
511 
512  if (okcol && okrow)
513  {
514  return kTRUE;
515  }
516  else
517  {
518  AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
519 
520  return kFALSE;
521  }
522 }
523 
534 //_______________________________________________________________________________
536  const UShort_t* cellList,
537  Int_t nCells)
538 {
539  if (!fRemoveBadChannels) return kFALSE;
540  if (!fEMCALBadChannelMap) return kFALSE;
541 
542  Int_t icol = -1;
543  Int_t irow = -1;
544  Int_t imod = -1;
545  for (Int_t iCell = 0; iCell<nCells; iCell++)
546  {
547  //Get the column and row
548  Int_t iTower = -1, iIphi = -1, iIeta = -1;
549  geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
550 
551  if (fEMCALBadChannelMap->GetEntries() <= imod) continue;
552 
553  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
554 
555  if (GetEMCALChannelStatus(imod, icol, irow))
556  {
557  AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
558  return kTRUE;
559  }
560  }// cell cluster loop
561 
562  return kFALSE;
563 }
564 
576 //___________________________________________________________________________
578  AliVCaloCells* cells, Int_t bc)
579 {
580  AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
581 
582  if(!geom)
583  {
584  AliError("No instance of the geometry is available");
585  return -1;
586  }
587 
588  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
589  geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
590  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
591 
592  // Get close cells index, energy and time, not in corners
593 
594  Int_t absID1 = -1;
595  Int_t absID2 = -1;
596 
597  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
598  if ( iphi > 0 ) absID2 = geom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
599 
600  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
601 
602  Int_t absID3 = -1;
603  Int_t absID4 = -1;
604 
605  if ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) )
606  {
607  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
608  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
609  }
610  else if ( ieta == 0 && imod%2 )
611  {
612  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
613  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
614  }
615  else
616  {
617  if ( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
618  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
619  if ( ieta > 0 )
620  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
621  }
622 
623  //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4);
624 
625  Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
626  Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
627 
628  AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
629  AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
630  AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
631  AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
632 
633  if (TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
634  if (TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
635  if (TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
636  if (TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
637 
638  return ecell1+ecell2+ecell3+ecell4;
639 }
640 
651 //_____________________________________________________________________________________________
652 Bool_t AliEMCALRecoUtils::IsExoticCell(Int_t absID, AliVCaloCells* cells, Int_t bc)
653 {
654  if (!fRejectExoticCells) return kFALSE;
655 
656  Float_t ecell = 0;
657  Double_t tcell = 0;
658  Bool_t accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
659 
660  if (!accept) return kTRUE; // reject this cell
661 
662  if (ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
663 
664  Float_t eCross = GetECross(absID,tcell,cells,bc);
665 
666  if (1-eCross/ecell > fExoticCellFraction)
667  {
668  AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",
669  absID,ecell,eCross,1-eCross/ecell));
670  return kTRUE;
671  }
672 
673  return kFALSE;
674 }
675 
685 //___________________________________________________________________
686 Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster,
687  AliVCaloCells *cells,
688  Int_t bc)
689 {
690  if (!cluster)
691  {
692  AliInfo("Cluster pointer null!");
693  return kFALSE;
694  }
695 
696  if (!fRejectExoticCluster) return kFALSE;
697 
698  // Get highest energy tower
699  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
700 
701  if(!geom)
702  {
703  AliError("No instance of the geometry is available");
704  return kFALSE;
705  }
706 
707  Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
708  Bool_t shared = kFALSE;
709  GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
710 
711  return IsExoticCell(absId,cells,bc);
712 }
713 
722 //_______________________________________________________________________
723 Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster)
724 {
725  if (!cluster)
726  {
727  AliInfo("Cluster pointer null!");
728  return 0;
729  }
730 
731  Float_t energy = cluster->E() ;
732  Float_t rdmEnergy = energy ;
733 
734  if (fSmearClusterEnergy)
735  {
736  rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
737  fSmearClusterParam[1] * energy +
738  fSmearClusterParam[2] );
739  AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
740  }
741 
742  return rdmEnergy;
743 }
744 
753 //____________________________________________________________________________
755 {
756  if (!cluster)
757  {
758  AliInfo("Cluster pointer null!");
759  return 0;
760  }
761 
762  Float_t energy = cluster->E();
763  if (energy < 0.100)
764  {
765  // Clusters with less than 100 MeV or negative are not possible
766  AliInfo(Form("Too low cluster energy!, E = %f < 0.100 GeV",energy));
767  return 0;
768  }
769  else if(energy > 300.)
770  {
771  AliInfo(Form("Too high cluster energy!, E = %f GeV, do not apply non linearity",energy));
772  return energy;
773  }
774 
775  switch (fNonLinearityFunction)
776  {
777  case kPi0MC:
778  {
779  //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
780  //fNonLinearityParams[0] = 1.014;
781  //fNonLinearityParams[1] =-0.03329;
782  //fNonLinearityParams[2] =-0.3853;
783  //fNonLinearityParams[3] = 0.5423;
784  //fNonLinearityParams[4] =-0.4335;
785  energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
786  ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
787  exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
788  break;
789  }
790 
791  case kPi0MCv2:
792  {
793  //Non-Linearity correction (from MC with function [0]/((x+[1])^[2]))+1;
794  //fNonLinearityParams[0] = 3.11111e-02;
795  //fNonLinearityParams[1] =-5.71666e-02;
796  //fNonLinearityParams[2] = 5.67995e-01;
797 
798  energy *= fNonLinearityParams[0]/TMath::Power(energy+fNonLinearityParams[1],fNonLinearityParams[2])+1;
799  break;
800  }
801 
802  case kPi0MCv3:
803  {
804  //Same as beam test corrected, change parameters
805  //fNonLinearityParams[0] = 9.81039e-01
806  //fNonLinearityParams[1] = 1.13508e-01;
807  //fNonLinearityParams[2] = 1.00173e+00;
808  //fNonLinearityParams[3] = 9.67998e-02;
809  //fNonLinearityParams[4] = 2.19381e+02;
810  //fNonLinearityParams[5] = 6.31604e+01;
811  //fNonLinearityParams[6] = 1;
812  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
813 
814  break;
815  }
816 
817 
818  case kPi0GammaGamma:
819  {
820  //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
821  //fNonLinearityParams[0] = 1.04;
822  //fNonLinearityParams[1] = -0.1445;
823  //fNonLinearityParams[2] = 1.046;
824  energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
825  break;
826  }
827 
828  case kPi0GammaConversion:
829  {
830  //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
831  //fNonLinearityParams[0] = 0.139393/0.1349766;
832  //fNonLinearityParams[1] = 0.0566186;
833  //fNonLinearityParams[2] = 0.982133;
834  energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
835 
836  break;
837  }
838 
839  case kBeamTest:
840  {
841  //From beam test, Alexei's results, for different ZS thresholds
842  // th=30 MeV; th = 45 MeV; th = 75 MeV
843  //fNonLinearityParams[0] = 1.007; 1.003; 1.002
844  //fNonLinearityParams[1] = 0.894; 0.719; 0.797
845  //fNonLinearityParams[2] = 0.246; 0.334; 0.358
846  //Rescale the param[0] with 1.03
847  energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
848 
849  break;
850  }
851 
852  case kBeamTestCorrected:
853  {
854  // From beam test, corrected for material between beam and EMCAL
855  //fNonLinearityParams[0] = 0.99078
856  //fNonLinearityParams[1] = 0.161499;
857  //fNonLinearityParams[2] = 0.655166;
858  //fNonLinearityParams[3] = 0.134101;
859  //fNonLinearityParams[4] = 163.282;
860  //fNonLinearityParams[5] = 23.6904;
861  //fNonLinearityParams[6] = 0.978;
862  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
863 
864  break;
865  }
866 
868  {
869  // From beam test, corrected for material between beam and EMCAL
870  // Different function to kBeamTestCorrected
871  //fNonLinearityParams[0] = 0.983504;
872  //fNonLinearityParams[1] = 0.210106;
873  //fNonLinearityParams[2] = 0.897274;
874  //fNonLinearityParams[3] = 0.0829064;
875  //fNonLinearityParams[4] = 152.299;
876  //fNonLinearityParams[5] = 31.5028;
877  //fNonLinearityParams[6] = 0.968;
878  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
879 
880  break;
881  }
882 
884  {
885  // Same function as kBeamTestCorrectedv2, different default parametrization.
886  //fNonLinearityParams[0] = 0.976941;
887  //fNonLinearityParams[1] = 0.162310;
888  //fNonLinearityParams[2] = 1.08689;
889  //fNonLinearityParams[3] = 0.0819592;
890  //fNonLinearityParams[4] = 152.338;
891  //fNonLinearityParams[5] = 30.9594;
892  //fNonLinearityParams[6] = 0.9615;
893  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
894 
895  break;
896  }
897 
898  case kSDMv5:
899  {
900  //Based on fit to the MC/data using kNoCorrection on the data - utilizes symmetric decay method and kPi0MCv5(MC) - 28 Oct 2013
901  //fNonLinearityParams[0] = 1.0;
902  //fNonLinearityParams[1] = 6.64778e-02;
903  //fNonLinearityParams[2] = 1.570;
904  //fNonLinearityParams[3] = 9.67998e-02;
905  //fNonLinearityParams[4] = 2.19381e+02;
906  //fNonLinearityParams[5] = 6.31604e+01;
907  //fNonLinearityParams[6] = 1.01286;
908  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5])))) * (0.964 + exp(-3.132-0.435*energy*2.0));
909 
910  break;
911  }
912 
913  case kPi0MCv5:
914  {
915  //Based on comparing MC truth information to the reconstructed energy of clusters.
916  //fNonLinearityParams[0] = 1.0;
917  //fNonLinearityParams[1] = 6.64778e-02;
918  //fNonLinearityParams[2] = 1.570;
919  //fNonLinearityParams[3] = 9.67998e-02;
920  //fNonLinearityParams[4] = 2.19381e+02;
921  //fNonLinearityParams[5] = 6.31604e+01;
922  //fNonLinearityParams[6] = 1.01286;
923  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
924 
925  break;
926  }
927 
928  case kSDMv6:
929  {
930  //Based on fit to the MC/data using kNoCorrection on the data
931  // - utilizes symmetric decay method and kPi0MCv6(MC) - 09 Dec 2014
932  // - parameters constrained by the test beam data as well
933  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
934  //fNonLinearityParams[0] = 1.0;
935  //fNonLinearityParams[1] = 0.237767;
936  //fNonLinearityParams[2] = 0.651203;
937  //fNonLinearityParams[3] = 0.183741;
938  //fNonLinearityParams[4] = 155.427;
939  //fNonLinearityParams[5] = 17.0335;
940  //fNonLinearityParams[6] = 0.987054;
941  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
942 
943  break;
944  }
945 
946  case kPi0MCv6:
947  {
948  //Based on comparing MC truth information to the reconstructed energy of clusters.
949  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
950  //fNonLinearityParams[0] = 1.0;
951  //fNonLinearityParams[1] = 0.0797873;
952  //fNonLinearityParams[2] = 1.68322;
953  //fNonLinearityParams[3] = 0.0806098;
954  //fNonLinearityParams[4] = 244.586;
955  //fNonLinearityParams[5] = 116.938;
956  //fNonLinearityParams[6] = 1.00437;
957  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
958 
959  break;
960  }
961 
962  case kPCMv1:
963  {
964  //based on symmetric decays of pi0 meson
965  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
966  // parameters vary from MC to MC
967  //fNonLinearityParams[0] = 0.984876;
968  //fNonLinearityParams[1] = -9.999609;
969  //fNonLinearityParams[2] = -4.999891;
970  //fNonLinearityParams[3] = 0.;
971  //fNonLinearityParams[4] = 0.;
972  //fNonLinearityParams[5] = 0.;
973  //fNonLinearityParams[6] = 0.;
974  energy /= TMath::Power(fNonLinearityParams[0] + exp(fNonLinearityParams[1] + fNonLinearityParams[2]*energy),2);//result coming from calo-conv method
975 
976  break;
977  }
978 
979  case kPCMplusBTCv1:
980  {
981  //convolution of TestBeamCorrectedv3 with PCM method
982  //Based on comparing MC truth information to the reconstructed energy of clusters.
983  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
984  // parameters vary from MC to MC
985  //fNonLinearityParams[0] = 0.976941;
986  //fNonLinearityParams[1] = 0.162310;
987  //fNonLinearityParams[2] = 1.08689;
988  //fNonLinearityParams[3] = 0.0819592;
989  //fNonLinearityParams[4] = 152.338;
990  //fNonLinearityParams[5] = 30.9594;
991  //fNonLinearityParams[6] = 0.9615;
992  //fNonLinearityParams[7] = 0.984876;
993  //fNonLinearityParams[8] = -9.999609;
994  //fNonLinearityParams[9] = -4.999891;
995  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
996  energy /= TMath::Power(fNonLinearityParams[7] + exp(fNonLinearityParams[8] + fNonLinearityParams[9]*energy),2);//result coming from calo-conv method
997 
998  break;
999  }
1000 
1001  case kPCMsysv1:
1002  {
1003  // Systematic variation of kPCMv1
1004  //Based on comparing MC truth information to the reconstructed energy of clusters.
1005  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1006  // parameters vary from MC to MC
1007  //fNonLinearityParams[0] = 0.0;
1008  //fNonLinearityParams[1] = 1.0;
1009  //fNonLinearityParams[2] = 1.0;
1010  //fNonLinearityParams[3] = 0.0;
1011  //fNonLinearityParams[4] = 1.0;
1012  //fNonLinearityParams[5] = 0.0;
1013  //fNonLinearityParams[6] = 0.0;
1014  energy /= TMath::Power( (fNonLinearityParams[0] + fNonLinearityParams[1] * TMath::Power(energy,fNonLinearityParams[2]) ) /
1015  (fNonLinearityParams[3] + fNonLinearityParams[4] * TMath::Power(energy,fNonLinearityParams[5]) ) + fNonLinearityParams[6], 2);//result coming from calo-conv method
1016 
1017  break;
1018  }
1019 
1020 
1021  case kNoCorrection:
1022  AliDebug(2,"No correction on the energy\n");
1023  break;
1024 
1025  }
1026 
1027  return energy;
1028 }
1029 
1034 //__________________________________________________
1036 {
1037  if (fNonLinearityFunction == kPi0MC) {
1038  fNonLinearityParams[0] = 1.014;
1039  fNonLinearityParams[1] = -0.03329;
1040  fNonLinearityParams[2] = -0.3853;
1041  fNonLinearityParams[3] = 0.5423;
1042  fNonLinearityParams[4] = -0.4335;
1043  }
1044 
1046  fNonLinearityParams[0] = 3.11111e-02;
1047  fNonLinearityParams[1] =-5.71666e-02;
1048  fNonLinearityParams[2] = 5.67995e-01;
1049  }
1050 
1052  fNonLinearityParams[0] = 9.81039e-01;
1053  fNonLinearityParams[1] = 1.13508e-01;
1054  fNonLinearityParams[2] = 1.00173e+00;
1055  fNonLinearityParams[3] = 9.67998e-02;
1056  fNonLinearityParams[4] = 2.19381e+02;
1057  fNonLinearityParams[5] = 6.31604e+01;
1058  fNonLinearityParams[6] = 1;
1059  }
1060 
1062  fNonLinearityParams[0] = 1.04;
1063  fNonLinearityParams[1] = -0.1445;
1064  fNonLinearityParams[2] = 1.046;
1065  }
1066 
1068  fNonLinearityParams[0] = 0.139393;
1069  fNonLinearityParams[1] = 0.0566186;
1070  fNonLinearityParams[2] = 0.982133;
1071  }
1072 
1074  if (fNonLinearThreshold == 30) {
1075  fNonLinearityParams[0] = 1.007;
1076  fNonLinearityParams[1] = 0.894;
1077  fNonLinearityParams[2] = 0.246;
1078  }
1079  if (fNonLinearThreshold == 45) {
1080  fNonLinearityParams[0] = 1.003;
1081  fNonLinearityParams[1] = 0.719;
1082  fNonLinearityParams[2] = 0.334;
1083  }
1084  if (fNonLinearThreshold == 75) {
1085  fNonLinearityParams[0] = 1.002;
1086  fNonLinearityParams[1] = 0.797;
1087  fNonLinearityParams[2] = 0.358;
1088  }
1089  }
1090 
1092  fNonLinearityParams[0] = 0.99078;
1093  fNonLinearityParams[1] = 0.161499;
1094  fNonLinearityParams[2] = 0.655166;
1095  fNonLinearityParams[3] = 0.134101;
1096  fNonLinearityParams[4] = 163.282;
1097  fNonLinearityParams[5] = 23.6904;
1098  fNonLinearityParams[6] = 0.978;
1099  }
1100 
1102  // Parameters until November 2015, use now kBeamTestCorrectedv3
1103  fNonLinearityParams[0] = 0.983504;
1104  fNonLinearityParams[1] = 0.210106;
1105  fNonLinearityParams[2] = 0.897274;
1106  fNonLinearityParams[3] = 0.0829064;
1107  fNonLinearityParams[4] = 152.299;
1108  fNonLinearityParams[5] = 31.5028;
1109  fNonLinearityParams[6] = 0.968;
1110  }
1111 
1113 
1114  // New parametrization of kBeamTestCorrectedv2
1115  // excluding point at 0.5 GeV from Beam Test Data
1116  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1117 
1118  fNonLinearityParams[0] = 0.976941;
1119  fNonLinearityParams[1] = 0.162310;
1120  fNonLinearityParams[2] = 1.08689;
1121  fNonLinearityParams[3] = 0.0819592;
1122  fNonLinearityParams[4] = 152.338;
1123  fNonLinearityParams[5] = 30.9594;
1124  fNonLinearityParams[6] = 0.9615;
1125  }
1126 
1127  if (fNonLinearityFunction == kSDMv5) {
1128  fNonLinearityParams[0] = 1.0;
1129  fNonLinearityParams[1] = 6.64778e-02;
1130  fNonLinearityParams[2] = 1.570;
1131  fNonLinearityParams[3] = 9.67998e-02;
1132  fNonLinearityParams[4] = 2.19381e+02;
1133  fNonLinearityParams[5] = 6.31604e+01;
1134  fNonLinearityParams[6] = 1.01286;
1135  }
1136 
1138  fNonLinearityParams[0] = 1.0;
1139  fNonLinearityParams[1] = 6.64778e-02;
1140  fNonLinearityParams[2] = 1.570;
1141  fNonLinearityParams[3] = 9.67998e-02;
1142  fNonLinearityParams[4] = 2.19381e+02;
1143  fNonLinearityParams[5] = 6.31604e+01;
1144  fNonLinearityParams[6] = 1.01286;
1145  }
1146 
1147  if (fNonLinearityFunction == kSDMv6) {
1148  fNonLinearityParams[0] = 1.0;
1149  fNonLinearityParams[1] = 0.237767;
1150  fNonLinearityParams[2] = 0.651203;
1151  fNonLinearityParams[3] = 0.183741;
1152  fNonLinearityParams[4] = 155.427;
1153  fNonLinearityParams[5] = 17.0335;
1154  fNonLinearityParams[6] = 0.987054;
1155  }
1156 
1158  fNonLinearityParams[0] = 1.0;
1159  fNonLinearityParams[1] = 0.0797873;
1160  fNonLinearityParams[2] = 1.68322;
1161  fNonLinearityParams[3] = 0.0806098;
1162  fNonLinearityParams[4] = 244.586;
1163  fNonLinearityParams[5] = 116.938;
1164  fNonLinearityParams[6] = 1.00437;
1165  }
1166 
1167 if (fNonLinearityFunction == kPCMv1) {
1168  //parameters change from MC production to MC production, they need to set for each period
1169  fNonLinearityParams[0] = 0.984876;
1170  fNonLinearityParams[1] = -9.999609;
1171  fNonLinearityParams[2] = -4.999891;
1172  fNonLinearityParams[3] = 0.;
1173  fNonLinearityParams[4] = 0.;
1174  fNonLinearityParams[5] = 0.;
1175  fNonLinearityParams[6] = 0.;
1176  }
1177 
1179  // test beam corrected values convoluted with symmetric meson decays values
1180  // for test beam:
1181  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1182  // for PCM method:
1183  // https://aliceinfo.cern.ch/Notes/node/211
1184  fNonLinearityParams[0] = 0.976941;
1185  fNonLinearityParams[1] = 0.162310;
1186  fNonLinearityParams[2] = 1.08689;
1187  fNonLinearityParams[3] = 0.0819592;
1188  fNonLinearityParams[4] = 152.338;
1189  fNonLinearityParams[5] = 30.9594;
1190  fNonLinearityParams[6] = 0.9615;
1191  fNonLinearityParams[7] = 0.984876;
1192  fNonLinearityParams[8] = -9.999609;
1193  fNonLinearityParams[9] = -4.999891;
1194  }
1195 
1197  //systematics for kPCMv1
1198  // for PCM method:
1199  // https://aliceinfo.cern.ch/Notes/node/211
1200  fNonLinearityParams[0] = 0.0;
1201  fNonLinearityParams[1] = 1.0;
1202  fNonLinearityParams[2] = 1.0;
1203  fNonLinearityParams[3] = 0.0;
1204  fNonLinearityParams[4] = 1.0;
1205  fNonLinearityParams[5] = 0.0;
1206  fNonLinearityParams[6] = 0.0;
1207  }
1208 }
1209 
1210 
1224 //____________________________________________________________________
1225 void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
1226  AliVCaloCells* cells,
1227  const AliVCluster* clu,
1228  Int_t & absId,
1229  Int_t & iSupMod,
1230  Int_t & ieta,
1231  Int_t & iphi,
1232  Bool_t & shared)
1233 {
1234  Double_t eMax = -1.;
1235  Double_t eCell = -1.;
1236  Float_t fraction = 1.;
1237  Float_t recalFactor = 1.;
1238  Int_t cellAbsId = -1 ;
1239 
1240  Int_t iTower = -1;
1241  Int_t iIphi = -1;
1242  Int_t iIeta = -1;
1243  Int_t iSupMod0= -1;
1244 
1245  if (!clu)
1246  {
1247  AliInfo("Cluster pointer null!");
1248  absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
1249  return;
1250  }
1251 
1252  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1253  {
1254  cellAbsId = clu->GetCellAbsId(iDig);
1255  fraction = clu->GetCellAmplitudeFraction(iDig);
1256  //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
1257  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1258 
1259  geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1260  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1261 
1262  if (iDig==0)
1263  {
1264  iSupMod0=iSupMod;
1265  } else if (iSupMod0!=iSupMod)
1266  {
1267  shared = kTRUE;
1268  //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
1269  }
1270 
1272  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1273 
1274  eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
1275  //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
1276  if (eCell > eMax)
1277  {
1278  eMax = eCell;
1279  absId = cellAbsId;
1280  //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
1281  }
1282  }// cell loop
1283 
1284  //Get from the absid the supermodule, tower and eta/phi numbers
1285  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1286  //Gives SuperModule and Tower numbers
1287  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1288  //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
1289  //printf("Max end---\n");
1290 }
1291 
1300 //_________________________________________________________
1302 {
1303  if (eCell > 0 && eCluster > 0)
1304  {
1305  if ( fW0 > 0 ) return TMath::Max( 0., fW0 + TMath::Log( eCell / eCluster ) ) ;
1306  else return TMath::Log( eCluster / eCell ) ;
1307  }
1308  else
1309  return 0. ;
1310 }
1311 
1315 //______________________________________
1317 {
1318  fParticleType = kPhoton;
1319  fPosAlgo = kUnchanged;
1320  fW0 = 4.5;
1321 
1323  fNonLinearThreshold = 30;
1324 
1325  fExoticCellFraction = 0.97;
1326  fExoticCellDiffTime = 1e6;
1328 
1329  fAODFilterMask = 128;
1330  fAODHybridTracks = kFALSE;
1331  fAODTPCOnlyTracks = kTRUE;
1332 
1333  fCutEtaPhiSum = kTRUE;
1334  fCutEtaPhiSeparate = kFALSE;
1335 
1336  fCutR = 0.05;
1337  fCutEta = 0.025;
1338  fCutPhi = 0.05;
1339 
1340  fClusterWindow = 100;
1341  fMass = 0.139;
1342 
1343  fStepSurface = 20.;
1344  fStepCluster = 5.;
1346 
1347  fCutMinTrackPt = 0;
1348  fCutMinNClusterTPC = -1;
1349  fCutMinNClusterITS = -1;
1350 
1351  fCutMaxChi2PerClusterTPC = 1e10;
1352  fCutMaxChi2PerClusterITS = 1e10;
1353 
1354  fCutRequireTPCRefit = kFALSE;
1355  fCutRequireITSRefit = kFALSE;
1356  fCutAcceptKinkDaughters = kFALSE;
1357 
1358  fCutMaxDCAToVertexXY = 1e10;
1359  fCutMaxDCAToVertexZ = 1e10;
1360  fCutDCAToVertex2D = kFALSE;
1361 
1362  fCutRequireITSStandAlone = kFALSE; //MARCEL
1363  fCutRequireITSpureSA = kFALSE; //Marcel
1364 
1365  // Misalignment matrices
1366  for (Int_t i = 0; i < 15 ; i++)
1367  {
1368  fMisalTransShift[i] = 0.;
1369  fMisalRotShift[i] = 0.;
1370  }
1371 
1372  // Non linearity
1373  for (Int_t i = 0; i < 10 ; i++) fNonLinearityParams[i] = 0.;
1374 
1375  // For kBeamTestCorrectedv2 case, but default is no correction
1376  fNonLinearityParams[0] = 0.983504;
1377  fNonLinearityParams[1] = 0.210106;
1378  fNonLinearityParams[2] = 0.897274;
1379  fNonLinearityParams[3] = 0.0829064;
1380  fNonLinearityParams[4] = 152.299;
1381  fNonLinearityParams[5] = 31.5028;
1382  fNonLinearityParams[6] = 0.968;
1383 
1384  // Cluster energy smearing
1385  fSmearClusterEnergy = kFALSE;
1386  fSmearClusterParam[0] = 0.07; // * sqrt E term
1387  fSmearClusterParam[1] = 0.00; // * E term
1388  fSmearClusterParam[2] = 0.00; // constant
1389 }
1390 
1394 //_____________________________________________________
1396 {
1397  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1398 
1399  // In order to avoid rewriting the same histograms
1400  Bool_t oldStatus = TH1::AddDirectoryStatus();
1401  TH1::AddDirectory(kFALSE);
1402 
1404  for (int i = 0; i < 22; i++)
1405  fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1406  Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
1407  //Init the histograms with 1
1408  for (Int_t sm = 0; sm < 22; sm++)
1409  {
1410  for (Int_t i = 0; i < 48; i++)
1411  {
1412  for (Int_t j = 0; j < 24; j++)
1413  {
1415  }
1416  }
1417  }
1418 
1419  fEMCALRecalibrationFactors->SetOwner(kTRUE);
1420  fEMCALRecalibrationFactors->Compress();
1421 
1422  // In order to avoid rewriting the same histograms
1423  TH1::AddDirectory(oldStatus);
1424 }
1425 
1429 //_________________________________________________________
1431 {
1432  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1433 
1434  // In order to avoid rewriting the same histograms
1435  Bool_t oldStatus = TH1::AddDirectoryStatus();
1436  TH1::AddDirectory(kFALSE);
1437 
1440 
1441  for (int i = 0; i < 4; i++)
1442  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1443  Form("hAllTimeAvBC%d",i),
1444  48*24*22,0.,48*24*22) );
1445  // Init the histograms with 0
1446  for (Int_t iBC = 0; iBC < 4; iBC++)
1447  {
1448  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1449  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kFALSE);
1450  }
1451 
1452  if(fLowGain) {
1453  for (int iBC = 0; iBC < 4; iBC++) {
1454  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvLGBC%d",iBC),
1455  Form("hAllTimeAvLGBC%d",iBC),
1456  48*24*22,0.,48*24*22) );
1457  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1458  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kTRUE);
1459  }
1460  }
1461 
1462 
1463 
1464  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1465  fEMCALTimeRecalibrationFactors->Compress();
1466 
1467  // In order to avoid rewriting the same histograms
1468  TH1::AddDirectory(oldStatus);
1469 }
1470 
1474 //____________________________________________________
1476 {
1477  AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1478 
1479  // In order to avoid rewriting the same histograms
1480  Bool_t oldStatus = TH1::AddDirectoryStatus();
1481  TH1::AddDirectory(kFALSE);
1482 
1483  fEMCALBadChannelMap = new TObjArray(22);
1484  //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
1485 
1486  for (int i = 0; i < 22; i++)
1487  fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1488 
1489  fEMCALBadChannelMap->SetOwner(kTRUE);
1490  fEMCALBadChannelMap->Compress();
1491 
1492  // In order to avoid rewriting the same histograms
1493  TH1::AddDirectory(oldStatus);
1494 }
1495 
1499 //___________________________________________________________
1501 {
1502  AliDebug(2,"AliEMCALRecoUtils::InitEMCALL1PhaseInTimeRecalibrationFactors()");
1503 
1504  // In order to avoid rewriting the same histograms
1505  Bool_t oldStatus = TH1::AddDirectoryStatus();
1506  TH1::AddDirectory(kFALSE);
1507 
1509 
1510  fEMCALL1PhaseInTimeRecalibration->Add(new TH1C("h0","EMCALL1phaseForSM", 22, 0, 22));
1511 
1512  for (Int_t i = 0; i < 22; i++) //loop over SMs, default value = 0
1514 
1515  fEMCALL1PhaseInTimeRecalibration->SetOwner(kTRUE);
1517 
1518  // In order to avoid rewriting the same histograms
1519  TH1::AddDirectory(oldStatus);
1520 }
1521 
1531 //____________________________________________________________________________
1532 void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
1533  AliVCluster * cluster,
1534  AliVCaloCells * cells,
1535  Int_t bc)
1536 {
1537  if (!cluster)
1538  {
1539  AliInfo("Cluster pointer null!");
1540  return;
1541  }
1542 
1543  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1544  UShort_t * index = cluster->GetCellsAbsId() ;
1545  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1546  Int_t ncells = cluster->GetNCells();
1547 
1548  // Initialize some used variables
1549  Float_t energy = 0;
1550  Int_t absId =-1;
1551  Int_t icol =-1, irow =-1, imod=1;
1552  Float_t factor = 1, frac = 0;
1553  Int_t absIdMax = -1;
1554  Float_t emax = 0;
1555 
1556  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1557  for (Int_t icell = 0; icell < ncells; icell++)
1558  {
1559  absId = index[icell];
1560  frac = fraction[icell];
1561  if (frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
1562 
1564  {
1565  // Energy
1566  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1567  geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1568  if (fEMCALRecalibrationFactors->GetEntries() <= imod)
1569  continue;
1570  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1571  factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1572 
1573  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1574  imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1575 
1576  }
1577 
1578  energy += cells->GetCellAmplitude(absId)*factor*frac;
1579 
1580  if (emax < cells->GetCellAmplitude(absId)*factor*frac)
1581  {
1582  emax = cells->GetCellAmplitude(absId)*factor*frac;
1583  absIdMax = absId;
1584  }
1585  }
1586 
1587  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
1588 
1589  cluster->SetE(energy);
1590 
1591  // Recalculate time of cluster
1592  Double_t timeorg = cluster->GetTOF();
1593  Bool_t isLowGain = !(cells->GetCellHighGain(absIdMax));//HG = false -> LG = true
1594 
1595  Double_t time = cells->GetCellTime(absIdMax);
1597  RecalibrateCellTime(absIdMax,bc,time,isLowGain);
1598  time-=fConstantTimeShift*1e-9; // only in case of Run1 old simulations
1599 
1600  // Recalibrate time with L1 phase
1602  RecalibrateCellTimeL1Phase(imod, bc, time);
1603 
1604  cluster->SetTOF(time);
1605 
1606  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
1607 }
1608 
1616 //_______________________________________________________________________
1617 void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, Int_t bc)
1618 {
1620  return;
1621 
1622  if (!cells)
1623  {
1624  AliInfo("Cells pointer null!");
1625  return;
1626  }
1627 
1628  Short_t absId =-1;
1629  Bool_t accept = kFALSE;
1630  Float_t ecell = 0;
1631  Double_t tcell = 0;
1632  Double_t ecellin = 0;
1633  Double_t tcellin = 0;
1634  Int_t mclabel = -1;
1635  Double_t efrac = 0;
1636 
1637  Int_t nEMcell = cells->GetNumberOfCells() ;
1638  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1639  {
1640  cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
1641 
1642  accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
1643  if (!accept)
1644  {
1645  ecell = 0;
1646  tcell = -1;
1647  }
1648 
1649  // Set new values
1650  cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
1651  }
1652 
1653  fCellsRecalibrated = kTRUE;
1654 }
1655 
1664 //_______________________________________________________________________________________________________
1665 void AliEMCALRecoUtils::RecalibrateCellTime(Int_t absId, Int_t bc, Double_t & celltime, Bool_t isLGon) const
1666 {
1667  if (!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0) {
1668  if(fLowGain)
1669  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,isLGon)*1.e-9;
1670  else
1671  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,kFALSE)*1.e-9;
1672  }
1673 }
1674 
1682 //_______________________________________________________________________________________________________
1684 {
1685  if (!fCellsRecalibrated && IsL1PhaseInTimeRecalibrationOn() && bc >= 0)
1686  {
1687  bc=bc%4;
1688 
1689  Float_t offsetPerSM=0.;
1690  Int_t l1PhaseShift = GetEMCALL1PhaseInTimeRecalibrationForSM(iSM);
1691  Int_t l1Phase=l1PhaseShift & 3; //bit operation
1692 
1693  if(bc >= l1Phase)
1694  offsetPerSM = (bc - l1Phase)*25;
1695  else
1696  offsetPerSM = (bc - l1Phase + 4)*25;
1697 
1698  Int_t l1shiftOffset=l1PhaseShift>>2; //bit operation
1699  l1shiftOffset*=25;
1700 
1701  celltime -= offsetPerSM*1.e-9;
1702  celltime -= l1shiftOffset*1.e-9;
1703  }
1704 }
1705 
1721 //______________________________________________________________________________
1722 void AliEMCALRecoUtils::RecalculateCellLabelsRemoveAddedGenerator( Int_t absID, AliVCluster* clus, AliMCEvent* mc,
1723  Float_t & amp, TArrayI & labeArr, TArrayF & eDepArr) const
1724 {
1725  TString genName ;
1726  Float_t eDepFrac[4];
1727 
1728  Float_t edepTotFrac = 1;
1729  Bool_t found = kFALSE;
1730  Float_t ampOrg = amp;
1731 
1732  //
1733  // Get the energy deposition fraction from cluster.
1734  //
1735  for(Int_t icluscell = 0; icluscell < clus->GetNCells(); icluscell++ )
1736  {
1737  if(absID == clus->GetCellAbsId(icluscell))
1738  {
1739  clus->GetCellMCEdepFractionArray(icluscell,eDepFrac);
1740 
1741  found = kTRUE;
1742 
1743  break;
1744  }
1745  }
1746 
1747  if ( !found )
1748  {
1749  AliWarning(Form("Cell abs ID %d NOT found in cluster",absID));
1750  return;
1751  }
1752 
1753  //
1754  // Check if there is a particle contribution from a given generator name.
1755  // If it is not one of the selected generators,
1756  // remove the constribution from the cell.
1757  //
1758  if ( mc && fNMCGenerToAccept > 0 )
1759  {
1760  //printf("Accept contribution from generator? \n");
1761  for(UInt_t imc = 0; imc < 4; imc++)
1762  {
1763  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc )
1764  {
1765  mc->GetCocktailGenerator(clus->GetLabelAt(imc),genName);
1766 
1767  Bool_t generOK = kFALSE;
1768  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
1769  {
1770  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
1771  }
1772 
1773  if ( !generOK )
1774  {
1775  amp-=ampOrg*eDepFrac[imc];
1776 
1777  edepTotFrac-=eDepFrac[imc];
1778 
1779  eDepFrac[imc] = 0;
1780  }
1781 
1782  } // eDep > 0
1783  } // 4 possible loop
1784 
1785  } // accept at least one kind of generator
1786 
1787  //
1788  // Add MC label and Edep to corresponding array (to be used later in digits)
1789  //
1790  Int_t nLabels = 0;
1791  for(UInt_t imc = 0; imc < 4; imc++)
1792  {
1793  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc && edepTotFrac > 0 )
1794  {
1795  nLabels++;
1796 
1797  labeArr.Set(nLabels);
1798  labeArr.AddAt(clus->GetLabelAt(imc), nLabels-1);
1799 
1800  eDepArr.Set(nLabels);
1801  eDepArr.AddAt( (eDepFrac[imc]/edepTotFrac) * amp, nLabels-1);
1802  // use as deposited energy a fraction of the simulated energy (smeared and with noise)
1803  } // edep > 0
1804 
1805  } // mc cell label loop
1806 
1807  //
1808  // If no label found, reject cell
1809  // It can happen to have this case (4 MC labels per cell is not enough for some cases)
1810  // better to remove. To be treated carefully.
1811  //
1812  if ( nLabels == 0 ) amp = 0;
1813 
1814 }
1815 
1823 //______________________________________________________________________________
1824 void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
1825  AliVCaloCells* cells,
1826  AliVCluster* clu)
1827 {
1828  if (!clu)
1829  {
1830  AliInfo("Cluster pointer null!");
1831  return;
1832  }
1833 
1835  else if (fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1836  else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
1837 }
1838 
1847 //_____________________________________________________________________________________________
1849  AliVCaloCells* cells,
1850  AliVCluster* clu)
1851 {
1852  Double_t eCell = 0.;
1853  Float_t fraction = 1.;
1854  Float_t recalFactor = 1.;
1855 
1856  Int_t absId = -1;
1857  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1858  Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
1859  Float_t weight = 0., totalWeight=0.;
1860  Float_t newPos[3] = {-1.,-1.,-1.};
1861  Double_t pLocal[3], pGlobal[3];
1862  Bool_t shared = kFALSE;
1863 
1864  Float_t clEnergy = clu->E(); //Energy already recalibrated previously
1865  if (clEnergy <= 0) return;
1866 
1867  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
1868  Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1869 
1870  //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1871 
1872  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1873  {
1874  absId = clu->GetCellAbsId(iDig);
1875  fraction = clu->GetCellAmplitudeFraction(iDig);
1876  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1877 
1878  if (!fCellsRecalibrated)
1879  {
1880  geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
1881  geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
1882  if (IsRecalibrationOn()) {
1883  recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1884  }
1885  }
1886 
1887  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1888 
1889  weight = GetCellWeight(eCell,clEnergy);
1890  totalWeight += weight;
1891 
1892  geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
1893  //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
1894  geom->GetGlobal(pLocal,pGlobal,iSupModMax);
1895  //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1896 
1897  for (int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
1898  }// cell loop
1899 
1900  if (totalWeight>0)
1901  {
1902  for (int i=0; i<3; i++ ) newPos[i] /= totalWeight;
1903  }
1904 
1905  //Float_t pos[]={0,0,0};
1906  //clu->GetPosition(pos);
1907  //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
1908  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1909 
1910  if (iSupModMax > 1) { //sector 1
1911  newPos[0] +=fMisalTransShift[3];//-=3.093;
1912  newPos[1] +=fMisalTransShift[4];//+=6.82;
1913  newPos[2] +=fMisalTransShift[5];//+=1.635;
1914  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
1915  } else { //sector 0
1916  newPos[0] +=fMisalTransShift[0];//+=1.134;
1917  newPos[1] +=fMisalTransShift[1];//+=8.2;
1918  newPos[2] +=fMisalTransShift[2];//+=1.197;
1919  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
1920  }
1921  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1922 
1923  clu->SetPosition(newPos);
1924 }
1925 
1934 //____________________________________________________________________________________________
1936  AliVCaloCells* cells,
1937  AliVCluster* clu)
1938 {
1939  Double_t eCell = 1.;
1940  Float_t fraction = 1.;
1941  Float_t recalFactor = 1.;
1942 
1943  Int_t absId = -1;
1944  Int_t iTower = -1;
1945  Int_t iIphi = -1, iIeta = -1;
1946  Int_t iSupMod = -1, iSupModMax = -1;
1947  Int_t iphi = -1, ieta =-1;
1948  Bool_t shared = kFALSE;
1949 
1950  Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
1951 
1952  if (clEnergy <= 0)
1953  return;
1954  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
1955  Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1956 
1957  Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
1958  Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1959  Int_t startingSM = -1;
1960 
1961  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1962  {
1963  absId = clu->GetCellAbsId(iDig);
1964  fraction = clu->GetCellAmplitudeFraction(iDig);
1965  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1966 
1967  if (iDig==0) startingSM = iSupMod;
1968  else if (iSupMod != startingSM) areInSameSM = kFALSE;
1969 
1970  eCell = cells->GetCellAmplitude(absId);
1971 
1972  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1973  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1974 
1975  if (!fCellsRecalibrated)
1976  {
1977  if (IsRecalibrationOn()) {
1978  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1979  }
1980  }
1981 
1982  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1983 
1984  weight = GetCellWeight(eCell,clEnergy);
1985  if (weight < 0) weight = 0;
1986  totalWeight += weight;
1987  weightedCol += ieta*weight;
1988  weightedRow += iphi*weight;
1989 
1990  //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
1991  }// cell loop
1992 
1993  Float_t xyzNew[]={-1.,-1.,-1.};
1994  if (areInSameSM == kTRUE)
1995  {
1996  //printf("In Same SM\n");
1997  weightedCol = weightedCol/totalWeight;
1998  weightedRow = weightedRow/totalWeight;
1999  geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2000  }
2001  else
2002  {
2003  //printf("In Different SM\n");
2004  geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2005  }
2006 
2007  clu->SetPosition(xyzNew);
2008 }
2009 
2017 //___________________________________________________________________________________________
2019  AliVCaloCells* cells,
2020  AliVCluster * cluster)
2021 {
2022  if (!fRecalDistToBadChannels) return;
2023 
2024  if (!cluster)
2025  {
2026  AliInfo("Cluster pointer null!");
2027  return;
2028  }
2029 
2030  // Get channels map of the supermodule where the cluster is.
2031  Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
2032  Bool_t shared = kFALSE;
2033  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
2034  TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
2035 
2036  Int_t dRrow, dRcol;
2037  Float_t minDist = 10000.;
2038  Float_t dist = 0.;
2039 
2040  // Loop on tower status map
2041  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2042  {
2043  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2044  {
2045  // Check if tower is bad.
2046  if (hMap->GetBinContent(icol,irow)==0) continue;
2047  //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
2048  // iSupMod,icol, irow, icolM,irowM);
2049 
2050  dRrow=TMath::Abs(irowM-irow);
2051  dRcol=TMath::Abs(icolM-icol);
2052  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2053  if (dist < minDist)
2054  {
2055  //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
2056  minDist = dist;
2057  }
2058  }
2059  }
2060 
2061  // In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
2062  if (shared)
2063  {
2064  TH2D* hMap2 = 0;
2065  Int_t iSupMod2 = -1;
2066 
2067  // The only possible combinations are (0,1), (2,3) ... (8,9)
2068  if (iSupMod%2) iSupMod2 = iSupMod-1;
2069  else iSupMod2 = iSupMod+1;
2070  hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
2071 
2072  // Loop on tower status map of second super module
2073  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2074  {
2075  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2076  {
2077  // Check if tower is bad.
2078  if (hMap2->GetBinContent(icol,irow)==0) continue;
2079 
2080  //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",
2081  // iSupMod2,icol, irow,iSupMod,icolM,irowM);
2082  dRrow=TMath::Abs(irow-irowM);
2083 
2084  if (iSupMod%2)
2085  dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
2086  else
2087  dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
2088 
2089  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2090  if (dist < minDist) minDist = dist;
2091  }
2092  }
2093  }// shared cluster in 2 SuperModules
2094 
2095  AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
2096  cluster->SetDistanceToBadChannel(minDist);
2097 }
2098 
2105 //__________________________________________________________________
2106 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
2107 {
2108  if (!cluster)
2109  {
2110  AliInfo("Cluster pointer null!");
2111  return;
2112  }
2113 
2114  if (cluster->GetM02() != 0)
2115  fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
2116 
2117  Float_t pidlist[AliPID::kSPECIESCN+1];
2118  for (Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
2119 
2120  cluster->SetPID(pidlist);
2121 }
2122 
2143 //___________________________________________________________________________________________________________________
2145  AliVCaloCells* cells, AliVCluster * cluster,
2146  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2147  Float_t & enAfterCuts, Float_t & l0, Float_t & l1,
2148  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2149  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2150 {
2151  if (!cluster)
2152  {
2153  AliInfo("Cluster pointer null!");
2154  return;
2155  }
2156 
2157  Double_t eCell = 0.;
2158  Double_t tCell = 0.;
2159  Float_t fraction = 1.;
2160  Float_t recalFactor = 1.;
2161  Bool_t isLowGain = kFALSE;
2162 
2163  Int_t iSupMod = -1;
2164  Int_t iTower = -1;
2165  Int_t iIphi = -1;
2166  Int_t iIeta = -1;
2167  Int_t iphi = -1;
2168  Int_t ieta = -1;
2169  Double_t etai = -1.;
2170  Double_t phii = -1.;
2171 
2172  Int_t nstat = 0 ;
2173  Float_t wtot = 0.;
2174  Double_t w = 0.;
2175  Double_t etaMean = 0.;
2176  Double_t phiMean = 0.;
2177 
2178  Double_t pGlobal[3];
2179 
2180  // Loop on cells, calculate the cluster energy, in case a cut on cell energy is added,
2181  // or the non linearity correction was applied
2182  // and to check if the cluster is between 2 SM in eta
2183  Int_t iSM0 = -1;
2184  Bool_t shared = kFALSE;
2185  Float_t energy = 0;
2186 
2187  enAfterCuts = 0;
2188  l0 = 0; l1 = 0;
2189  disp = 0; dEta = 0; dPhi = 0;
2190  sEta = 0; sPhi = 0; sEtaPhi = 0;
2191 
2192  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2193  {
2194  // Get from the absid the supermodule, tower and eta/phi numbers
2195 
2196  Int_t absId = cluster->GetCellAbsId(iDigit);
2197 
2198  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2199  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2200 
2201  // Check if there are cells of different SM
2202  if (iDigit == 0 ) iSM0 = iSupMod;
2203  else if (iSupMod!= iSM0) shared = kTRUE;
2204 
2205  // Get the cell energy, if recalibration is on, apply factors
2206  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2207  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2208 
2209  if (IsRecalibrationOn())
2210  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2211 
2212  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2213  tCell = cells->GetCellTime (absId);
2214  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2215 
2216  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2217  tCell*=1e9;
2218  tCell-=fConstantTimeShift;
2219 
2220  if(eCell > 0.05) // at least a minimum 50 MeV cut
2221  energy += eCell;
2222 
2223  if(eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut)
2224  enAfterCuts += eCell;
2225  } // cell loop
2226 
2227 
2228  // Loop on cells to calculate weights and shower shape terms parameters
2229  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2230  {
2231  // Get from the absid the supermodule, tower and eta/phi numbers
2232  Int_t absId = cluster->GetCellAbsId(iDigit);
2233 
2234  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2235  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2236 
2237  //Get the cell energy, if recalibration is on, apply factors
2238  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2239  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2240 
2242  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2243 
2244  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2245  tCell = cells->GetCellTime (absId);
2246  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2247 
2248  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2249  tCell*=1e9;
2250  tCell-=fConstantTimeShift;
2251 
2252  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2253  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2254  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2255 
2256  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2257  {
2258  w = GetCellWeight(eCell, energy);
2259 
2260  // Cell index
2261  if ( fShowerShapeCellLocationType == 0 )
2262  {
2263  etai=(Double_t)ieta;
2264  phii=(Double_t)iphi;
2265  }
2266  // Cell angle location
2267  else if( fShowerShapeCellLocationType == 1 )
2268  {
2269  geom->EtaPhiFromIndex(absId, etai, phii);
2270  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2271  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2272  }
2273  else
2274  {
2275  geom->GetGlobal(absId,pGlobal);
2276 
2277  // Cell x-z location
2278  if( fShowerShapeCellLocationType == 2 )
2279  {
2280  etai = pGlobal[2];
2281  phii = pGlobal[0];
2282  }
2283  // Cell r-z location
2284  else
2285  {
2286  etai = pGlobal[2];
2287  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2288  }
2289  }
2290 
2291  if (w > 0.0)
2292  {
2293  wtot += w ;
2294  nstat++;
2295 
2296  // Shower shape
2297  sEta += w * etai * etai ;
2298  etaMean += w * etai ;
2299  sPhi += w * phii * phii ;
2300  phiMean += w * phii ;
2301  sEtaPhi += w * etai * phii ;
2302  }
2303  }
2304  else if(eCell > 0.05)
2305  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2306  } // cell loop
2307 
2308  //printf("sEta %f sPhi %f etaMean %f phiMean %f sEtaPhi %f wtot %f\n",sEta,sPhi,etaMean,phiMean,sEtaPhi, wtot);
2309 
2310  // Normalize to the weight
2311  if (wtot > 0)
2312  {
2313  etaMean /= wtot ;
2314  phiMean /= wtot ;
2315  }
2316  else
2317  AliDebug(2,Form("Wrong weight %f\n", wtot));
2318 
2319  // Loop on cells to calculate dispersion
2320  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2321  {
2322  // Get from the absid the supermodule, tower and eta/phi numbers
2323  Int_t absId = cluster->GetCellAbsId(iDigit);
2324 
2325  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2326  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2327 
2328  //Get the cell energy, if recalibration is on, apply factors
2329  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2330  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2331  if (IsRecalibrationOn())
2332  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2333 
2334  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2335  tCell = cells->GetCellTime (absId);
2336  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2337 
2338  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2339  tCell*=1e9;
2340  tCell-=fConstantTimeShift;
2341 
2342  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2343  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2344  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2345 
2346  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2347  {
2348  w = GetCellWeight(eCell,cluster->E());
2349 
2350  // Cell index
2351  if ( fShowerShapeCellLocationType == 0 )
2352  {
2353  etai=(Double_t)ieta;
2354  phii=(Double_t)iphi;
2355  }
2356  // Cell angle location
2357  else if( fShowerShapeCellLocationType == 1 )
2358  {
2359  geom->EtaPhiFromIndex(absId, etai, phii);
2360  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2361  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2362  }
2363  else
2364  {
2365  geom->GetGlobal(absId,pGlobal);
2366 
2367  // Cell x-z location
2368  if( fShowerShapeCellLocationType == 2 )
2369  {
2370  etai = pGlobal[2];
2371  phii = pGlobal[0];
2372  }
2373  // Cell r-z location
2374  else
2375  {
2376  etai = pGlobal[2];
2377  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2378  }
2379  }
2380 
2381  if (w > 0.0)
2382  {
2383  disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
2384  dEta += w * (etai-etaMean)*(etai-etaMean) ;
2385  dPhi += w * (phii-phiMean)*(phii-phiMean) ;
2386  }
2387  }
2388  else
2389  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2390  }// cell loop
2391 
2392  // Normalize to the weigth and set shower shape parameters
2393  if (wtot > 0 && nstat > 1)
2394  {
2395  disp /= wtot ;
2396  dEta /= wtot ;
2397  dPhi /= wtot ;
2398  sEta /= wtot ;
2399  sPhi /= wtot ;
2400  sEtaPhi /= wtot ;
2401 
2402  sEta -= etaMean * etaMean ;
2403  sPhi -= phiMean * phiMean ;
2404  sEtaPhi -= etaMean * phiMean ;
2405 
2406  l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2407  l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2408 
2409  //printf("sEta %f sPhi %f etaMean %f phiMean %f sEtaPhi %f wtot %f l0 %f l1 %f\n",sEta,sPhi,etaMean,phiMean,sEtaPhi, wtot,l0,l1);
2410 
2411  }
2412  else
2413  {
2414  l0 = 0. ;
2415  l1 = 0. ;
2416  dEta = 0. ; dPhi = 0. ; disp = 0. ;
2417  sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
2418  }
2419 }
2420 
2439 //___________________________________________________________________________________________________________________
2441  AliVCaloCells* cells, AliVCluster * cluster,
2442  Float_t & l0, Float_t & l1,
2443  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2444  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2445 {
2446  Float_t newEnergy = 0;
2447  Float_t cellEmin = 0.05; // 50 MeV
2448  Float_t cellTimeWindow = 1000; // open cut
2449  Int_t bc = 0;
2451  cellEmin, cellTimeWindow, bc,
2452  newEnergy, l0, l1, disp,
2453  dEta, dPhi, sEta, sPhi, sEtaPhi);
2454 }
2455 
2464 //____________________________________________________________________________________________
2466  AliVCaloCells* cells,
2467  AliVCluster * cluster)
2468 {
2469  Float_t l0 = 0., l1 = 0.;
2470  Float_t disp = 0., dEta = 0., dPhi = 0.;
2471  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2472 
2473  AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
2474  dEta, dPhi, sEta, sPhi, sEtaPhi);
2475 
2476  cluster->SetM02(l0);
2477  cluster->SetM20(l1);
2478  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2479 
2480 }
2481 
2494 //____________________________________________________________________________________________
2496  AliVCaloCells* cells, AliVCluster * cluster,
2497  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2498  Float_t & enAfterCuts)
2499 {
2500  Float_t l0 = 0., l1 = 0.;
2501  Float_t disp = 0., dEta = 0., dPhi = 0.;
2502  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2503 
2505  cellEcut, cellTimeCut, bc,
2506  enAfterCuts, l0, l1, disp,
2507  dEta, dPhi, sEta, sPhi, sEtaPhi);
2508 
2509  cluster->SetM02(l0);
2510  cluster->SetM20(l1);
2511  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2512 
2513 }
2514 
2528 //____________________________________________________________________________
2529 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
2530  TObjArray * clusterArr,
2531  const AliEMCALGeometry *geom,
2532  AliMCEvent * mc)
2533 {
2534  fMatchedTrackIndex ->Reset();
2535  fMatchedClusterIndex->Reset();
2536  fResidualPhi->Reset();
2537  fResidualEta->Reset();
2538 
2539  fMatchedTrackIndex ->Set(1000);
2540  fMatchedClusterIndex->Set(1000);
2541  fResidualPhi->Set(1000);
2542  fResidualEta->Set(1000);
2543 
2544  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
2545  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
2546 
2547  // Init the magnetic field if not already on
2548  if (!TGeoGlobalMagField::Instance()->GetField())
2549  {
2550  if (!event->InitMagneticField())
2551  {
2552  AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
2553  }
2554  } // Init mag field
2555 
2556  if (esdevent)
2557  {
2558  UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
2559  UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
2560  Bool_t desc1 = (mask1 >> 3) & 0x1;
2561  Bool_t desc2 = (mask2 >> 3) & 0x1;
2562  if (desc1==0 || desc2==0) {
2563  // AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
2564  // mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
2565  // mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
2566  fITSTrackSA=kTRUE;
2567  }
2568  }
2569 
2570  TObjArray *clusterArray = 0x0;
2571  if (!clusterArr)
2572  {
2573  clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
2574  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2575  {
2576  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2577  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells()))
2578  continue;
2579  clusterArray->AddAt(cluster,icl);
2580  }
2581  }
2582 
2583  Int_t matched=0;
2584  Double_t cv[21];
2585  TString genName;
2586  for (Int_t i=0; i<21;i++) cv[i]=0;
2587  for (Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
2588  {
2589  AliExternalTrackParam *trackParam = 0;
2590  Int_t mcLabel = -1;
2591  // If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
2592  AliESDtrack *esdTrack = 0;
2593  AliAODTrack *aodTrack = 0;
2594  if (esdevent)
2595  {
2596  esdTrack = esdevent->GetTrack(itr);
2597  if (!esdTrack) continue;
2598  if (!IsAccepted(esdTrack)) continue;
2599  if (esdTrack->Pt()<fCutMinTrackPt) continue;
2600 
2601  if ( TMath::Abs(esdTrack->Eta()) > 0.9 ) continue;
2602 
2603  // Save some time and memory in case of no DCal present
2604  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2605  {
2606  Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
2607  if ( phi <= 10 || phi >= 250 ) continue;
2608  }
2609 
2610  if (!fITSTrackSA)
2611  trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); // if TPC Available
2612  else
2613  trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
2614 
2615  mcLabel = TMath::Abs(esdTrack->GetLabel());
2616  }
2617 
2618  // If the input event is AOD, the starting point for extrapolation is at vertex
2619  // AOD tracks are selected according to its filterbit.
2620  else if (aodevent)
2621  {
2622  aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(itr));
2623 
2624  if (!aodTrack) AliFatal("Not a standard AOD");
2625 
2626  if (!aodTrack) continue;
2627 
2628  if (fAODTPCOnlyTracks)
2629  { // Match with TPC only tracks, default from May 2013, before filter bit 32
2630  //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
2631  if (!aodTrack->IsTPCConstrained()) continue ;
2632  }
2633  else if (fAODHybridTracks)
2634  { // Match with hybrid tracks
2635  //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
2636  if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
2637  } else
2638  { // Match with tracks on a mask
2639  //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask));
2640  if (!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks
2641  }
2642 
2643  if (aodTrack->Pt() < fCutMinTrackPt) continue;
2644 
2645  if ( TMath::Abs(aodTrack->Eta()) > 0.9 ) continue;
2646 
2647  // Save some time and memory in case of no DCal present
2648  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2649  {
2650  Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
2651  if ( phi <= 10 || phi >= 250 ) continue;
2652  }
2653 
2654  Double_t pos[3],mom[3];
2655  aodTrack->GetXYZ(pos);
2656  aodTrack->GetPxPyPz(mom);
2657  AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",
2658  itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
2659 
2660  trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
2661 
2662  mcLabel = TMath::Abs(aodTrack->GetLabel());
2663  }
2664 
2665  //Return if the input data is not "AOD" or "ESD"
2666  else
2667  {
2668  AliWarning("Wrong input data type! Should be \"AOD\" or \"ESD\" ");
2669  if (clusterArray)
2670  {
2671  clusterArray->Clear();
2672  delete clusterArray;
2673  }
2674  return;
2675  }
2676 
2677  if (!trackParam) continue;
2678 
2679  //
2680  // Check if track comes from a particular MC generator, do not include it
2681  // if it is not a selected one
2682  //
2683  if( mc && fMCGenerToAcceptForTrack && fNMCGenerToAccept > 0 )
2684  {
2685  mc->GetCocktailGenerator(mcLabel,genName);
2686 
2687  Bool_t generOK = kFALSE;
2688  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
2689  {
2690  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
2691  }
2692 
2693  if ( !generOK ) continue;
2694  }
2695 
2696  // Extrapolate the track to EMCal surface, see AliEMCALRecoUtilsBase
2697  AliExternalTrackParam emcalParam(*trackParam);
2698  Float_t eta, phi, pt;
2699  if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2700  {
2701  if (aodevent && trackParam) delete trackParam;
2702  if (fITSTrackSA && trackParam) delete trackParam;
2703  continue;
2704  }
2705 
2706  if ( TMath::Abs(eta) > 0.75 )
2707  {
2708  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2709  continue;
2710  }
2711 
2712  // Save some time and memory in case of no DCal present
2713  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2714  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2715  {
2716  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2717  continue;
2718  }
2719 
2720  //Find matched clusters
2721  Int_t index = -1;
2722  Float_t dEta = -999, dPhi = -999;
2723  if (!clusterArr)
2724  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
2725  else
2726  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr , dEta, dPhi);
2727 
2728 
2729  if (index>-1)
2730  {
2731  fMatchedTrackIndex ->AddAt(itr,matched);
2732  fMatchedClusterIndex ->AddAt(index,matched);
2733  fResidualEta ->AddAt(dEta,matched);
2734  fResidualPhi ->AddAt(dPhi,matched);
2735  matched++;
2736  }
2737 
2738  if (aodevent && trackParam) delete trackParam;
2739  if (fITSTrackSA && trackParam) delete trackParam;
2740  }//track loop
2741 
2742  if (clusterArray)
2743  {
2744  clusterArray->Clear();
2745  delete clusterArray;
2746  }
2747 
2748  AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
2749 
2750  fMatchedTrackIndex ->Set(matched);
2751  fMatchedClusterIndex ->Set(matched);
2752  fResidualPhi ->Set(matched);
2753  fResidualEta ->Set(matched);
2754 }
2755 
2767 //________________________________________________________________________________
2769  const AliVEvent *event,
2770  const AliEMCALGeometry *geom,
2771  Float_t &dEta, Float_t &dPhi)
2772 {
2773  Int_t index = -1;
2774 
2775  if ( TMath::Abs(track->Eta()) > 0.9 ) return index;
2776 
2777  // Save some time and memory in case of no DCal present
2778  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2779  {
2780  Double_t phiV = track->Phi()*TMath::RadToDeg();
2781  if ( phiV <= 10 || phiV >= 250 ) return index;
2782  }
2783 
2784  AliExternalTrackParam *trackParam = 0;
2785  if (!fITSTrackSA)
2786  trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam()); // If TPC
2787  else
2788  trackParam = new AliExternalTrackParam(*track);
2789 
2790  if (!trackParam) return index;
2791  AliExternalTrackParam emcalParam(*trackParam);
2792 
2793  Float_t eta, phi, pt;
2794  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2795  {
2796  if (fITSTrackSA) delete trackParam;
2797  return index;
2798  }
2799 
2800  if ( TMath::Abs(eta) > 0.75 )
2801  {
2802  if (fITSTrackSA) delete trackParam;
2803  return index;
2804  }
2805 
2806  // Save some time and memory in case of no DCal present
2807  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2808  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2809  {
2810  if (fITSTrackSA) delete trackParam;
2811  return index;
2812  }
2813 
2814  TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
2815 
2816  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2817  {
2818  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2819  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
2820  clusterArr->AddAt(cluster,icl);
2821  }
2822 
2823  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
2824  clusterArr->Clear();
2825  delete clusterArr;
2826  if (fITSTrackSA) delete trackParam;
2827 
2828  return index;
2829 }
2830 
2841 //_______________________________________________________________________________________________
2842 Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
2843  AliExternalTrackParam *trkParam,
2844  const TObjArray * clusterArr,
2845  Float_t &dEta, Float_t &dPhi)
2846 {
2847  dEta=-999, dPhi=-999;
2848  Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
2849  Int_t index = -1;
2850  Float_t tmpEta=-999, tmpPhi=-999;
2851 
2852  Double_t exPos[3] = {0.,0.,0.};
2853  if (!emcalParam->GetXYZ(exPos)) return index;
2854 
2855  Float_t clsPos[3] = {0.,0.,0.};
2856  for (Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
2857  {
2858  AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
2859 
2860  if (!cluster || !cluster->IsEMCAL()) continue;
2861 
2862  cluster->GetPosition(clsPos);
2863 
2864  Double_t dR = TMath::Sqrt(TMath::Power(exPos[0]-clsPos[0],2)+TMath::Power(exPos[1]-clsPos[1],2)+TMath::Power(exPos[2]-clsPos[2],2));
2865  if (dR > fClusterWindow) continue;
2866 
2867  AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
2868 
2869  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
2870 
2871  if (fCutEtaPhiSum)
2872  {
2873  Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
2874  if (tmpR<dRMax)
2875  {
2876  dRMax=tmpR;
2877  dEtaMax=tmpEta;
2878  dPhiMax=tmpPhi;
2879  index=icl;
2880  }
2881  }
2882  else if (fCutEtaPhiSeparate)
2883  {
2884  if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
2885  {
2886  dEtaMax = tmpEta;
2887  dPhiMax = tmpPhi;
2888  index=icl;
2889  }
2890  }
2891  else
2892  {
2893  AliError("Please specify your cut criteria");
2894  AliError("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()");
2895  AliError("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()");
2896  return index;
2897  }
2898  }
2899 
2900  dEta=dEtaMax;
2901  dPhi=dPhiMax;
2902 
2903  return index;
2904 }
2905 
2918 //---------------------------------------------------------------------------------
2919 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
2920  const AliVCluster *cluster,
2921  Float_t &tmpEta,
2922  Float_t &tmpPhi)
2923 {
2924  return AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
2925 }
2926 
2936 //_______________________________________________________________________
2938  Float_t &dEta, Float_t &dPhi)
2939 {
2940  if (FindMatchedPosForCluster(clsIndex) >= 999)
2941  {
2942  AliDebug(2,"No matched tracks found!\n");
2943  dEta=999.;
2944  dPhi=999.;
2945  return;
2946  }
2947 
2948  dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
2949  dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
2950 }
2951 
2961 //______________________________________________________________________________________________
2963 {
2964  if (FindMatchedPosForTrack(trkIndex) >= 999)
2965  {
2966  AliDebug(2,"No matched cluster found!\n");
2967  dEta=999.;
2968  dPhi=999.;
2969  return;
2970  }
2971 
2972  dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
2973  dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
2974 }
2975 
2982 //__________________________________________________________
2984 {
2985  if (IsClusterMatched(clsIndex))
2986  return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
2987  else
2988  return -1;
2989 }
2990 
2997 //__________________________________________________________
2999 {
3000  if (IsTrackMatched(trkIndex))
3001  return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
3002  else
3003  return -1;
3004 }
3005 
3013 //______________________________________________________________
3015 {
3016  if (FindMatchedPosForCluster(clsIndex) < 999)
3017  return kTRUE;
3018  else
3019  return kFALSE;
3020 }
3021 
3029 //____________________________________________________________
3031 {
3032  if (FindMatchedPosForTrack(trkIndex) < 999)
3033  return kTRUE;
3034  else
3035  return kFALSE;
3036 }
3037 
3045 //______________________________________________________________________
3047 {
3048  Float_t tmpR = fCutR;
3049  UInt_t pos = 999;
3050 
3051  for (Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
3052  {
3053  if (fMatchedClusterIndex->At(i)==clsIndex)
3054  {
3055  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3056  if (r<tmpR)
3057  {
3058  pos=i;
3059  tmpR=r;
3060  AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3061  fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3062  }
3063  }
3064  }
3065  return pos;
3066 }
3067 
3075 //____________________________________________________________________
3077 {
3078  Float_t tmpR = fCutR;
3079  UInt_t pos = 999;
3080 
3081  for (Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
3082  {
3083  if (fMatchedTrackIndex->At(i)==trkIndex)
3084  {
3085  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3086  if (r<tmpR)
3087  {
3088  pos=i;
3089  tmpR=r;
3090  AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3091  fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3092  }
3093  }
3094  }
3095  return pos;
3096 }
3097 
3112 //__________________________________________________________________________
3114  const AliEMCALGeometry *geom,
3115  AliVCaloCells* cells, Int_t bc)
3116 {
3117  Bool_t isGood=kTRUE;
3118 
3119  if (!cluster || !cluster->IsEMCAL()) return kFALSE;
3120  if (ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
3121  if (!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
3122  if (IsExoticCluster(cluster, cells,bc)) return kFALSE;
3123 
3124  return isGood;
3125 }
3126 
3137 //__________________________________________________________
3139 {
3140  UInt_t status = esdTrack->GetStatus();
3141 
3142  Int_t nClustersITS = esdTrack->GetITSclusters(0);
3143  Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
3144 
3145  Float_t chi2PerClusterITS = -1;
3146  Float_t chi2PerClusterTPC = -1;
3147  if (nClustersITS!=0)
3148  chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
3149  if (nClustersTPC!=0)
3150  chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
3151 
3152  //
3153  // DCA cuts
3154  // Only to be used for primary
3155  //
3156  if ( fTrackCutsType == kGlobalCut )
3157  {
3158  Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01);
3159  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
3160 
3161  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3162 
3163  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3164  }
3165  else if( fTrackCutsType == kGlobalCut2011 )
3166  {
3167  Float_t maxDCAToVertexXYPtDep = 0.0105 + 0.0350/TMath::Power(esdTrack->Pt(),1.1);
3168  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2011()
3169 
3170  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3171 
3172  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3173  }
3174 
3175  Float_t b[2];
3176  Float_t bCov[3];
3177  esdTrack->GetImpactParameters(b,bCov);
3178  if (bCov[0]<=0 || bCov[2]<=0)
3179  {
3180  AliDebug(1, "Estimated b resolution lower or equal zero!");
3181  bCov[0]=0; bCov[2]=0;
3182  }
3183 
3184  Float_t dcaToVertexXY = b[0];
3185  Float_t dcaToVertexZ = b[1];
3186  Float_t dcaToVertex = -1;
3187 
3188  if (fCutDCAToVertex2D)
3189  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY / fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY +
3190  dcaToVertexZ*dcaToVertexZ / fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ );
3191  else
3192  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
3193 
3194  // cut the track?
3195 
3196  Bool_t cuts[kNCuts];
3197  for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
3198 
3199  // track quality cuts
3200  if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
3201  cuts[0]=kTRUE;
3202  if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
3203  cuts[1]=kTRUE;
3204  if (nClustersTPC<fCutMinNClusterTPC)
3205  cuts[2]=kTRUE;
3206  if (nClustersITS<fCutMinNClusterITS)
3207  cuts[3]=kTRUE;
3208  if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
3209  cuts[4]=kTRUE;
3210  if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
3211  cuts[5]=kTRUE;
3212  if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
3213  cuts[6]=kTRUE;
3214  if (fCutDCAToVertex2D && dcaToVertex > 1)
3215  cuts[7] = kTRUE;
3216  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
3217  cuts[8] = kTRUE;
3218  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
3219  cuts[9] = kTRUE;
3220 
3222  {
3223  //Require at least one SPD point + anything else in ITS
3224  if ( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
3225  cuts[10] = kTRUE;
3226  }
3227 
3228  // ITS
3230  {
3231  if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin))
3232  {
3233  // TPC tracks
3234  cuts[11] = kTRUE;
3235  }
3236  else
3237  {
3238  // ITS standalone tracks
3240  {
3241  if (status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
3242  }
3243  else if (fCutRequireITSpureSA)
3244  {
3245  if (!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
3246  }
3247  }
3248  }
3249 
3250  Bool_t cut=kFALSE;
3251  for (Int_t i=0; i<kNCuts; i++)
3252  if (cuts[i]) { cut = kTRUE ; }
3253 
3254  // cut the track
3255  if (cut)
3256  return kFALSE;
3257  else
3258  return kTRUE;
3259 }
3260 
3266 //_____________________________________
3268 {
3269  switch (fTrackCutsType)
3270  {
3271  case kTPCOnlyCut:
3272  {
3273  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()"));
3274  //TPC
3275  SetMinNClustersTPC(70);
3277  SetAcceptKinkDaughters(kFALSE);
3278  SetRequireTPCRefit(kFALSE);
3279 
3280  //ITS
3281  SetRequireITSRefit(kFALSE);
3282  SetMaxDCAToVertexZ(3.2);
3283  SetMaxDCAToVertexXY(2.4);
3284  SetDCAToVertex2D(kTRUE);
3285 
3286  break;
3287  }
3288 
3289  case kGlobalCut:
3290  {
3291  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE)"));
3292  //TPC
3293  SetMinNClustersTPC(70);
3295  SetAcceptKinkDaughters(kFALSE);
3296  SetRequireTPCRefit(kTRUE);
3297 
3298  //ITS
3299  SetRequireITSRefit(kTRUE);
3300  SetMaxDCAToVertexZ(2);
3302  SetDCAToVertex2D(kFALSE);
3303 
3304  break;
3305  }
3306 
3307  case kLooseCut:
3308  {
3309  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
3310  SetMinNClustersTPC(50);
3311  SetAcceptKinkDaughters(kTRUE);
3312 
3313  break;
3314  }
3315 
3316  case kITSStandAlone:
3317  {
3318  AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
3319  SetRequireITSRefit(kTRUE);
3320  SetRequireITSStandAlone(kTRUE);
3321  SetITSTrackSA(kTRUE);
3322  break;
3323  }
3324 
3325  case kGlobalCut2011:
3326  {
3327  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE)"));
3328  //TPC
3329  SetMinNClustersTPC(50);
3331  SetAcceptKinkDaughters(kFALSE);
3332  SetRequireTPCRefit(kTRUE);
3333 
3334  //ITS
3335  SetRequireITSRefit(kTRUE);
3336  SetMaxDCAToVertexZ(2);
3338  SetDCAToVertex2D(kFALSE);
3339 
3340  break;
3341  }
3342 
3343  case kLooseCutWithITSrefit:
3344  {
3345  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut plus ITSrefit"));
3346  SetMinNClustersTPC(50);
3347  SetAcceptKinkDaughters(kTRUE);
3348  SetRequireITSRefit(kTRUE);
3349 
3350  break;
3351  }
3352  }
3353 }
3354 
3360 //________________________________________________________________________
3362 {
3363  Int_t nTracks = event->GetNumberOfTracks();
3364  for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
3365  {
3366  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
3367  if (!track)
3368  {
3369  AliWarning(Form("Could not receive track %d", iTrack));
3370  continue;
3371  }
3372 
3373  Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
3374  track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
3375 
3376  /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
3377  AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
3378  if (esdtrack)
3379  {
3380  if (matchClusIndex != -1)
3381  esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
3382  else
3383  esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3384  }
3385  else
3386  {
3387  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
3388  if (matchClusIndex != -1)
3389  aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
3390  else
3391  aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3392  }
3393  }
3394  AliDebug(2,"Track matched to closest cluster");
3395 }
3396 
3403 //_________________________________________________________________________
3405 {
3406  for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
3407  {
3408  AliVCluster *cluster = event->GetCaloCluster(iClus);
3409  if (!cluster->IsEMCAL())
3410  continue;
3411 
3412  //
3413  // Remove old matches in cluster
3414  //
3415  if(cluster->GetNTracksMatched() > 0)
3416  {
3417  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
3418  {
3419  TArrayI arrayTrackMatched(0);
3420  ((AliESDCaloCluster*)cluster)->AddTracksMatched(arrayTrackMatched);
3421  }
3422  else
3423  {
3424  for(Int_t iTrack = 0; iTrack < cluster->GetNTracksMatched(); iTrack++)
3425  {
3426  ((AliAODCaloCluster*)cluster)->RemoveTrackMatched((TObject*)((AliAODCaloCluster*)cluster)->GetTrackMatched(iTrack));
3427  }
3428  }
3429  }
3430 
3431  //
3432  // Find new matches and put them in the cluster
3433  //
3434  Int_t nTracks = event->GetNumberOfTracks();
3435  TArrayI arrayTrackMatched(nTracks);
3436 
3437  // Get the closest track matched to the cluster
3438  Int_t nMatched = 0;
3439  Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
3440  if (matchTrackIndex != -1)
3441  {
3442  arrayTrackMatched[nMatched] = matchTrackIndex;
3443  nMatched++;
3444  }
3445 
3446  // Get all other tracks matched to the cluster
3447  for (Int_t iTrk=0; iTrk<nTracks; ++iTrk)
3448  {
3449  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
3450 
3451  if( !track ) continue;
3452 
3453  if ( iTrk == matchTrackIndex ) continue;
3454 
3455  if ( track->GetEMCALcluster() == iClus )
3456  {
3457  arrayTrackMatched[nMatched] = iTrk;
3458  ++nMatched;
3459  }
3460  }
3461 
3462  AliDebug(2,Form("cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]));
3463 
3464  arrayTrackMatched.Set(nMatched);
3465  AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
3466  if (esdcluster)
3467  esdcluster->AddTracksMatched(arrayTrackMatched);
3468  else if ( nMatched > 0 )
3469  {
3470  AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
3471  if (aodcluster)
3472  {
3473  aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
3474  //AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(event->GetTrack(arrayTrackMatched.At(0)));
3475  //printf("Is the closest matching track with ID %d a 128? %d what's its full filter map? %u\n",aodtrack->GetID(), aodtrack->TestFilterBit(128),aodtrack->GetFilterMap());
3476  //printf("With specs: pt %.4f, eta %.4f, phi %.4f\n",aodtrack->Pt(),aodtrack->Eta(), aodtrack->Phi());
3477  }
3478  }
3479 
3480  Float_t eta= -999, phi = -999;
3481  if (matchTrackIndex != -1)
3482  GetMatchedResiduals(iClus, eta, phi);
3483 
3484  cluster->SetTrackDistance(phi, eta);
3485  }
3486 
3487  AliDebug(2,"Cluster matched to tracks");
3488 }
3489 
3492  else {
3493  fEMCALRecalibrationFactors = new TObjArray(map->GetEntries());
3494  fEMCALRecalibrationFactors->SetOwner(true);
3495  }
3496  if(!fEMCALRecalibrationFactors->IsOwner()){
3497  // Must claim ownership since the new objects are owend by this instance
3498  fEMCALRecalibrationFactors->SetOwner(kTRUE);
3499  }
3500  for(int i = 0; i < map->GetEntries(); i++){
3501  TH2F *hist = dynamic_cast<TH2F *>(map->At(i));
3502  if(!hist) continue;
3503  this->SetEMCALChannelRecalibrationFactors(i, hist);
3504  }
3505 }
3506 
3510  fEMCALRecalibrationFactors->SetOwner(true);
3511  }
3512  if(fEMCALRecalibrationFactors->GetEntries() <= iSM) fEMCALRecalibrationFactors->Expand(iSM+1);
3513  if(fEMCALRecalibrationFactors->At(iSM)) fEMCALRecalibrationFactors->RemoveAt(iSM);
3514  TH2F *clone = new TH2F(*h);
3515  clone->SetDirectory(NULL);
3516  fEMCALRecalibrationFactors->AddAt(clone,iSM);
3517 }
3518 
3521  else {
3522  fEMCALBadChannelMap = new TObjArray(map->GetEntries());
3523  fEMCALBadChannelMap->SetOwner(true);
3524  }
3525  if(!fEMCALBadChannelMap->IsOwner()) {
3526  // Must claim ownership since the new objects are owend by this instance
3527  fEMCALBadChannelMap->SetOwner(true);
3528  }
3529  for(int i = 0; i < map->GetEntries(); i++){
3530  TH2I *hist = dynamic_cast<TH2I *>(map->At(i));
3531  if(!hist) continue;
3532  this->SetEMCALChannelStatusMap(i, hist);
3533  }
3534 }
3535 
3537  if(!fEMCALBadChannelMap){
3538  fEMCALBadChannelMap = new TObjArray(iSM);
3539  fEMCALBadChannelMap->SetOwner(true);
3540  }
3541  if(fEMCALBadChannelMap->GetEntries() <= iSM) fEMCALBadChannelMap->Expand(iSM+1);
3542  if(fEMCALBadChannelMap->At(iSM)) fEMCALBadChannelMap->RemoveAt(iSM);
3543  TH2I *clone = new TH2I(*h);
3544  clone->SetDirectory(NULL);
3545  fEMCALBadChannelMap->AddAt(clone,iSM);
3546 }
3547 
3550  else {
3551  fEMCALTimeRecalibrationFactors = new TObjArray(map->GetEntries());
3552  fEMCALTimeRecalibrationFactors->SetOwner(true);
3553  }
3554  if(!fEMCALTimeRecalibrationFactors->IsOwner()) {
3555  // Must claim ownership since the new objects are owend by this instance
3556  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
3557  }
3558  for(int i = 0; i < map->GetEntries(); i++){
3559  TH1F *hist = dynamic_cast<TH1F *>(map->At(i));
3560  if(!hist) continue;
3562  }
3563 }
3564 
3568  fEMCALTimeRecalibrationFactors->SetOwner(true);
3569  }
3570  if(fEMCALTimeRecalibrationFactors->GetEntries() <= bc) fEMCALTimeRecalibrationFactors->Expand(bc+1);
3572  TH1F *clone = new TH1F(*h);
3573  clone->SetDirectory(NULL);
3574  fEMCALTimeRecalibrationFactors->AddAt(clone,bc);
3575 }
3576 
3579  else {
3580  fEMCALL1PhaseInTimeRecalibration = new TObjArray(map->GetEntries());
3581  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3582  }
3583  if(!fEMCALL1PhaseInTimeRecalibration->IsOwner()) {
3584  // Must claim ownership since the new objects are owend by this instance
3585  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3586  }
3587  for(int i = 0; i < map->GetEntries(); i++) {
3588  TH1C *hist = dynamic_cast<TH1C *>(map->At(i));
3589  if(!hist) continue;
3591  }
3592 }
3593 
3597  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3598  }
3600  TH1C *clone = new TH1C(*h);
3601  clone->SetDirectory(NULL);
3602  fEMCALL1PhaseInTimeRecalibration->AddAt(clone,0);
3603 }
3604 
3608 //___________________________________________________
3610 {
3611  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3612  printf("AliEMCALRecoUtils Settings: \n");
3613  printf("\tMisalignment shifts\n");
3614  for (Int_t i=0; i<5; i++) printf("\t\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i,
3616  fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
3617  printf("\tNon linearity function %d, parameters:\n", fNonLinearityFunction);
3618  if (fNonLinearityFunction != 3) // print only if not kNoCorrection
3619  for (Int_t i=0; i<10; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
3620 
3621  printf("\tPosition Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
3622 
3623  printf("\tMatching criteria: ");
3624  if (fCutEtaPhiSum) {
3625  printf("\t\tsqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
3626  } else if (fCutEtaPhiSeparate) {
3627  printf("\t\tdEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
3628  } else {
3629  printf("\t\tError\n");
3630  printf("\t\tplease specify your cut criteria\n");
3631  printf("\t\tTo cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
3632  printf("\t\tTo cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
3633  }
3634 
3635  printf("\tMass hypothesis = %2.3f [GeV/c^2], extrapolation step to surface = %2.2f[cm], step to cluster = %2.2f[cm]\n",fMass,fStepSurface, fStepCluster);
3636  printf("\tCluster selection window: dR < %2.0f\n",fClusterWindow);
3637 
3638  printf("\tTrack cuts: \n");
3639  printf("\t\tMinimum track pT: %1.2f\n",fCutMinTrackPt);
3640  printf("\t\tAOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask);
3641  printf("\t\tTPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
3642  printf("\t\tAcceptKinks = %d\n",fCutAcceptKinkDaughters);
3643  printf("\t\tMinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
3644  printf("\t\tMaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
3645  printf("\t\tDCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
3646  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3647 }
void RecalculateCellLabelsRemoveAddedGenerator(Int_t absID, AliVCluster *clus, AliMCEvent *mc, Float_t &amp, TArrayI &labeArr, TArrayF &eDepArr) const
Bool_t IsL1PhaseInTimeRecalibrationOn() const
Double_t fMass
Mass hypothesis of the track.
Bool_t fMCGenerToAcceptForTrack
Activate the removal of tracks entering the track matching that come from a particular generator...
double Double_t
Definition: External.C:58
Bool_t CheckCellFiducialRegion(const AliEMCALGeometry *geom, const AliVCluster *cluster, AliVCaloCells *cells)
Bool_t fAODTPCOnlyTracks
Match with TPC only tracks.
Bool_t fCutRequireTPCRefit
Require TPC refit.
Definition: External.C:236
void SetRequireITSStandAlone(Bool_t b=kFALSE)
Bool_t IsClusterMatched(Int_t clsIndex) const
Float_t GetEMCALChannelTimeRecalibrationFactor(Int_t bc, Int_t absID, Bool_t isLGon=kFALSE) const
TArrayI * fMatchedClusterIndex
Array that stores indexes of matched clusters.
TObjArray * fEMCALTimeRecalibrationFactors
Array of histograms with map of time recalibration factors, EMCAL.
void SetMaxDCAToVertexXY(Float_t dist=1e10)
TRandom3 fRandom
Random generator for cluster energy smearing.
Bool_t fTimeRecalibration
Switch on or off the time recalibration.
Bool_t fCutEtaPhiSum
Place cut on sqrt(dEta^2+dPhi^2)
Int_t FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam, AliExternalTrackParam *trkParam, const TObjArray *clusterArr, Float_t &dEta, Float_t &dPhi)
void InitEMCALL1PhaseInTimeRecalibration()
Int_t fCutMinNClusterTPC
Min number of tpc clusters.
Float_t fCutMaxDCAToVertexZ
Track-to-vertex cut in max absolute distance in z-plane.
void SetRequireITSRefit(Bool_t b=kFALSE)
energy
Definition: HFPtSpectrum.C:44
void SetEMCALL1PhaseInTimeRecalibrationForSM(Int_t iSM, Int_t c=0)
void SetDCAToVertex2D(Bool_t b=kFALSE)
void SetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow, Double_t c=1)
Int_t fNonLinearThreshold
Non linearity threshold value for kBeamTest non linearity function.
Bool_t fITSTrackSA
If track matching is to be done with ITS tracks standing alone.
void SetITSTrackSA(Bool_t isITS)
Bool_t IsExoticCluster(const AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=0)
Float_t fMisalRotShift[15]
Cluster position rotation shift parameters.
void SetMinNClustersTPC(Int_t min=-1)
Bool_t fRecalDistToBadChannels
Calculate distance from highest energy tower of cluster to closes bad channel.
Bool_t fRecalibration
Switch on or off the recalibration.
TArrayI * fMatchedTrackIndex
Array that stores indexes of matched tracks.
void SetMaxDCAToVertexZ(Float_t dist=1e10)
Bool_t IsBadChannelsRemovalSwitchedOn() const
Bool_t IsExoticCell(Int_t absId, AliVCaloCells *cells, Int_t bc=-1)
void RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
Float_t fExoticCellDiffTime
If time of candidate to exotic and close cell is too different (in ns), it must be noisy...
TArrayF * fResidualPhi
Array that stores the residual phi.
void RecalibrateCells(AliVCaloCells *cells, Int_t bc)
Double_t fCutMinTrackPt
Cut on track pT.
Bool_t IsRecalibrationOn() const
Some utilities for cluster and cell treatment.
Float_t fCutMaxChi2PerClusterITS
Max its fit chi2 per its cluster.
void RecalculateClusterShowerShapeParametersWithCellCuts(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster, Float_t cellEcut, Float_t cellTimeCut, Int_t bc, Float_t &enAfterCuts)
Bool_t fRemoveBadChannels
Check the channel status provided and remove clusters with bad channels.
Bool_t IsAccepted(AliESDtrack *track)
Float_t fMisalTransShift[15]
Cluster position translation shift parameters.
void SetAcceptKinkDaughters(Bool_t b=kTRUE)
void RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
Double_t fClusterWindow
Select clusters in the window to be matched.
Bool_t fCutRequireITSRefit
Require ITS refit.
Float_t fNonLinearityParams[10]
Parameters for the non linearity function.
Float_t GetCellWeight(Float_t eCell, Float_t eCluster) const
Bool_t fUseL1PhaseInTimeRecalibration
Switch on or off the L1 phase in time recalibration.
Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow) const
void SetClusterMatchedToTrack(const AliVEvent *event)
Bool_t fCutDCAToVertex2D
If true a 2D DCA cut is made.
Int_t GetMatchedTrackIndex(Int_t clsIndex)
void RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
int Int_t
Definition: External.C:63
void Print(const Option_t *) const
UInt_t FindMatchedPosForCluster(Int_t clsIndex) const
Bool_t AcceptCalibrateCell(Int_t absId, Int_t bc, Float_t &amp, Double_t &time, AliVCaloCells *cells)
unsigned int UInt_t
Definition: External.C:33
void RecalculateClusterShowerShapeParameters(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
float Float_t
Definition: External.C:68
Float_t fCutR
sqrt(dEta^2+dPhi^2) cut on matching
Int_t fNCellsFromEMCALBorder
Number of cells from EMCAL border the cell with maximum amplitude has to be.
Bool_t fCellsRecalibrated
Internal bool to check if cells (time/energy) where recalibrated and not recalibrate them when recalc...
Bool_t fRejectExoticCluster
Switch on or off exotic cluster rejection.
Bool_t IsTimeRecalibrationOn() const
Int_t fParticleType
Particle type for depth calculation, see enum ParticleType.
Int_t GetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow) const
Definition: External.C:228
void SetMaxChi2PerClusterTPC(Float_t max=1e10)
Float_t fConstantTimeShift
Apply a 600 ns (+15.8) time shift in case of simulation, shift in ns.
Int_t GetMatchedClusterIndex(Int_t trkIndex)
Int_t fNonLinearityFunction
Non linearity function choice, see enum NonlinearityFunctions.
Bool_t fCutRequireITSStandAlone
Require ITSStandAlone.
TObjArray * fEMCALL1PhaseInTimeRecalibration
Histogram with map of L1 phase per SM, EMCAL.
TObjArray * fEMCALRecalibrationFactors
Array of histograms with map of recalibration factors, EMCAL.
Bool_t fUseRunCorrectionFactors
Use Run Dependent Correction.
void SetTracksMatchedToCluster(const AliVEvent *event)
Float_t fCutMaxChi2PerClusterTPC
Max tpc fit chi2 per tpc cluster.
Float_t fExoticCellMinAmplitude
Check for exotic only if amplitud is larger than this value.
void RecalibrateClusterEnergy(const AliEMCALGeometry *geom, AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=-1)
Float_t fW0
Energy weight used in cluster position and shower shape calculations.
Bool_t fAODHybridTracks
Match with hybrid.
void GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
AliEMCALRecoUtils & operator=(const AliEMCALRecoUtils &)
void SetRequireTPCRefit(Bool_t b=kFALSE)
TArrayF * fResidualEta
Array that stores the residual eta.
short Short_t
Definition: External.C:23
TString fMCGenerToAccept[5]
List with name of generators that should not be included.
Bool_t fSmearClusterEnergy
Smear cluster energy, to be done only for simulated data to match real data.
Bool_t IsGoodCluster(AliVCluster *cluster, const AliEMCALGeometry *geom, AliVCaloCells *cells, Int_t bc=-1)
void GetMaxEnergyCell(const AliEMCALGeometry *geom, AliVCaloCells *cells, const AliVCluster *clu, Int_t &absId, Int_t &iSupMod, Int_t &ieta, Int_t &iphi, Bool_t &shared)
void InitEMCALTimeRecalibrationFactors()
Double_t fStepSurface
Length of step to extrapolate tracks to EMCal surface.
UInt_t fAODFilterMask
Filter mask to select AOD tracks. Refer to $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C.
static Bool_t ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, const AliVCluster *cluster, Double_t mass, Double_t step, Float_t &tmpEta, Float_t &tmpPhi)
void GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
Float_t fSmearClusterParam[3]
Energy smearing parameters.
Bool_t fCutAcceptKinkDaughters
Accepting kink daughters?
void RecalculateClusterPosition(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
Float_t fExoticCellFraction
Good cell if fraction < 1-ecross/ecell.
void RecalculateClusterPID(AliVCluster *cluster)
Int_t fPosAlgo
Position recalculation algorithm, see enum PositionAlgorithms.
void SetEMCALL1PhaseInTimeRecalibrationForAllSM(const TObjArray *map)
AliEMCALPIDUtils * fPIDUtils
Recalculate PID parameters.
Bool_t fRejectExoticCells
Remove exotic cells.
Int_t FindMatchedClusterInEvent(const AliESDtrack *track, const AliVEvent *event, const AliEMCALGeometry *geom, Float_t &dEta, Float_t &dPhi)
Bool_t IsTrackMatched(Int_t trkIndex) const
void RecalibrateCellTimeL1Phase(Int_t iSM, Int_t bc, Double_t &time) const
Float_t fCutMaxDCAToVertexXY
Track-to-vertex cut in max absolute distance in xy-plane.
void FindMatches(AliVEvent *event, TObjArray *clusterArr=0x0, const AliEMCALGeometry *geom=0x0, AliMCEvent *mc=0x0)
unsigned short UShort_t
Definition: External.C:28
TObjArray * fEMCALBadChannelMap
Array of histograms with map of bad channels, EMCAL.
void SetEMCALChannelTimeRecalibrationFactors(const TObjArray *map)
Double_t fStepCluster
Length of step to extrapolate tracks to clusters.
Bool_t fNoEMCALBorderAtEta0
Do fiducial cut in EMCAL region eta = 0?
const char Option_t
Definition: External.C:48
Int_t fNMCGenerToAccept
Number of MC generators that should not be included in analysis.
void SetEMCALChannelTimeRecalibrationFactor(Int_t bc, Int_t absID, Double_t c=0, Bool_t isLGon=kFALSE)
Bool_t fCutRequireITSpureSA
ITS pure standalone tracks.
Int_t GetEMCALL1PhaseInTimeRecalibrationForSM(Int_t iSM) const
bool Bool_t
Definition: External.C:53
Int_t fCutMinNClusterITS
Min number of its clusters.
Float_t CorrectClusterEnergyLinearity(AliVCluster *clu)
UInt_t FindMatchedPosForTrack(Int_t trkIndex) const
Bool_t ClusterContainsBadChannel(const AliEMCALGeometry *geom, const UShort_t *cellList, Int_t nCells)
Double_t fEMCalSurfaceDistance
EMCal surface distance (= 430 by default, the last 10 cm are propagated on a cluster-track pair basis...
void SetEMCALChannelRecalibrationFactors(const TObjArray *map)
Float_t GetECross(Int_t absID, Double_t tcell, AliVCaloCells *cells, Int_t bc)
Bool_t fLowGain
Switch on or off calibration with low gain channels.
void SetEMCALChannelStatusMap(const TObjArray *map)
Float_t fCutPhi
dPhi cut on matching
Int_t fTrackCutsType
ESD track cuts type for matching, see enum TrackCutsType.
void RecalibrateCellTime(Int_t absId, Int_t bc, Double_t &time, Bool_t isLGon=kFALSE) const
Bool_t fCutEtaPhiSeparate
Cut on dEta and dPhi separately.
Float_t fCutEta
dEta cut on matching
Float_t SmearClusterEnergy(const AliVCluster *clu)