AliPhysics  master (3d17d9d)
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), fUseShaperNonlin(kFALSE),
56  fSmearClusterEnergy(kFALSE), fRandom(),
57  fCellsRecalibrated(kFALSE), fRecalibration(kFALSE), fUse1Drecalib(kFALSE), fEMCALRecalibrationFactors(),
58  fConstantTimeShift(0), fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(), fLowGain(kFALSE),
59  fUseL1PhaseInTimeRecalibration(kFALSE), fEMCALL1PhaseInTimeRecalibration(),
60  fIsParRun(kFALSE), fCurrentParNumber(0), fNPars(0), fGlobalEventID(NULL),
61  fDoUseMergedBC(kFALSE),
62  fUseRunCorrectionFactors(kFALSE),
63  fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(), fUse1Dmap(kFALSE),
64  fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
65  fRejectExoticCluster(kFALSE), fRejectExoticCells(kFALSE),
66  fExoticCellFraction(0), fExoticCellDiffTime(0), fExoticCellMinAmplitude(0),
67  fPIDUtils(), fAODFilterMask(0),
68  fAODHybridTracks(0), fAODTPCOnlyTracks(0),
69  fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
70  fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kFALSE), fCutEtaPhiSeparate(kFALSE),
71  fCutR(0), fCutEta(0), fCutPhi(0),
72  fClusterWindow(0), fMass(0),
73  fStepSurface(0), fStepCluster(0),
74  fITSTrackSA(kFALSE), fUseTrackDCA(kTRUE), // keep it active, but not working for old MC
75  fUseOuterTrackParam(kFALSE), fEMCalSurfaceDistance(440.),
76  fTrackCutsType(0), fCutMinTrackPt(0), fCutMinNClusterTPC(0),
77  fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
78  fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
79  fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0), fCutDCAToVertex2D(kFALSE),
80  fCutRequireITSStandAlone(kFALSE), fCutRequireITSpureSA(kFALSE),
81  fNMCGenerToAccept(0), fMCGenerToAcceptForTrack(1)
82 {
83  // Init parameters
85 
86  for(Int_t j = 0; j < 5; j++) fMCGenerToAccept[j] = "";
87 
88  // Track matching arrays init
91  fResidualPhi = new TArrayF();
92  fResidualEta = new TArrayF();
93  fPIDUtils = new AliEMCALPIDUtils();
94 
95  fBadStatusSelection[0] = kTRUE;
96  fBadStatusSelection[1] = kTRUE;
97  fBadStatusSelection[2] = kTRUE;
98  fBadStatusSelection[3] = kTRUE;
99 }
100 
101 //
102 // Copy constructor.
103 //
104 //______________________________________________________________________
106 : AliEMCALRecoUtilsBase(reco),
117  fLowGain(reco.fLowGain),
120  fIsParRun(reco.fIsParRun),
122  fNPars(reco.fNPars),
136  fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0),
137  fResidualPhi( reco.fResidualPhi? new TArrayF(*reco.fResidualPhi):0x0),
139  fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
152 {
153  for (Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ;
154  fMisalTransShift[i] = reco.fMisalTransShift[i] ; }
155  for (Int_t i = 0; i < 10 ; i++){ fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
156  for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
157  for (Int_t j = 0; j < 5 ; j++) { fMCGenerToAccept[j] = reco.fMCGenerToAccept[j] ; }
158  for (Int_t j = 0; j < 4 ; j++) { fBadStatusSelection[j] = reco.fBadStatusSelection[j] ; }
159 
160  if(reco.fEMCALBadChannelMap) {
161  // Copy constructor - not taking ownership over calibration histograms
162  fEMCALBadChannelMap = new TObjArray(reco.fEMCALBadChannelMap->GetEntries());
163  fEMCALBadChannelMap->SetOwner(false);
164  for(int ism = 0; ism < reco.fEMCALBadChannelMap->GetEntries(); ism++) fEMCALBadChannelMap->AddAt(reco.fEMCALBadChannelMap->At(ism), ism);
165  }
166 
167  if(reco.fEMCALRecalibrationFactors) {
168  // Copy constructor - not taking ownership over calibration histograms
170  fEMCALRecalibrationFactors->SetOwner(false);
171  for(int ism = 0; ism < reco.fEMCALRecalibrationFactors->GetEntries(); ism++) fEMCALRecalibrationFactors->AddAt(reco.fEMCALRecalibrationFactors->At(ism), ism);
172  }
173 
175  // Copy constructor - not taking ownership over calibration histograms
177  fEMCALTimeRecalibrationFactors->SetOwner(false);
178  for(int ism = 0; ism < reco.fEMCALTimeRecalibrationFactors->GetEntries(); ism++) fEMCALTimeRecalibrationFactors->AddAt(reco.fEMCALTimeRecalibrationFactors->At(ism), ism);
179  }
180 }
181 
185 //______________________________________________________________________
187 {
188  if (this == &reco)return *this;
189  ((TNamed *)this)->operator=(reco);
190 
191  for (Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ;
192  fMisalRotShift[i] = reco.fMisalRotShift[i] ; }
193  for (Int_t i = 0; i < 10 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
194  for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
195 
197  fPosAlgo = reco.fPosAlgo;
198  fW0 = reco.fW0;
200 
205 
209 
212  fLowGain = reco.fLowGain;
213 
216 
217  fIsParRun = reco.fIsParRun;
219  fNPars = reco.fNPars;
221 
223 
225 
228  fUse1Dmap = reco.fUse1Dmap;
229 
232 
238 
239  fPIDUtils = reco.fPIDUtils;
240 
244 
247  fCutR = reco.fCutR;
248  fCutEta = reco.fCutEta;
249  fCutPhi = reco.fCutPhi;
251  fMass = reco.fMass;
252  fStepSurface = reco.fStepSurface;
253  fStepCluster = reco.fStepCluster;
254  fITSTrackSA = reco.fITSTrackSA;
255  fUseTrackDCA = reco.fUseTrackDCA;
258 
273 
276  for (Int_t j = 0; j < 5 ; j++)
277  fMCGenerToAccept[j] = reco.fMCGenerToAccept[j];
278 
279  //
280  // Assign or copy construct the different TArrays
281  //
282  if (reco.fResidualEta)
283  {
284  if (fResidualEta)
285  *fResidualEta = *reco.fResidualEta;
286  else
287  fResidualEta = new TArrayF(*reco.fResidualEta);
288  }
289  else
290  {
291  if (fResidualEta) delete fResidualEta;
292  fResidualEta = 0;
293  }
294 
295  if (reco.fResidualPhi)
296  {
297  if (fResidualPhi)
298  *fResidualPhi = *reco.fResidualPhi;
299  else
300  fResidualPhi = new TArrayF(*reco.fResidualPhi);
301  }
302  else
303  {
304  if (fResidualPhi) delete fResidualPhi;
305  fResidualPhi = 0;
306  }
307 
308  if (reco.fMatchedTrackIndex)
309  {
310  if (fMatchedTrackIndex)
312  else
314  }
315  else
316  {
318  fMatchedTrackIndex = 0;
319  }
320 
321  if (reco.fMatchedClusterIndex)
322  {
325  else
327  }
328  else
329  {
332  }
333 
334  for (Int_t j = 0; j < 4 ; j++)
336 
338  if(reco.fEMCALBadChannelMap) {
339  // Copy constructor - not taking ownership over calibration histograms
340  fEMCALBadChannelMap = new TObjArray(reco.fEMCALBadChannelMap->GetEntries());
341  fEMCALBadChannelMap->SetOwner(false);
342  for(int ism = 0; ism < reco.fEMCALBadChannelMap->GetEntries(); ism++) fEMCALBadChannelMap->AddAt(reco.fEMCALBadChannelMap->At(ism), ism);
343  }
344 
346  if(reco.fEMCALRecalibrationFactors) {
347  // Copy constructor - not taking ownership over calibration histograms
349  fEMCALRecalibrationFactors->SetOwner(false);
350  for(int ism = 0; ism < reco.fEMCALRecalibrationFactors->GetEntries(); ism++) fEMCALRecalibrationFactors->AddAt(reco.fEMCALRecalibrationFactors->At(ism), ism);
351  }
352 
355  // Copy constructor - not taking ownership over calibration histograms
357  fEMCALTimeRecalibrationFactors->SetOwner(false);
358  for(int ism = 0; ism < reco.fEMCALTimeRecalibrationFactors->GetEntries(); ism++) fEMCALTimeRecalibrationFactors->AddAt(reco.fEMCALTimeRecalibrationFactors->At(ism), ism);
359  }
360 
363  // Copy constructor - not taking ownership over calibration histograms
365  fEMCALL1PhaseInTimeRecalibration->SetOwner(false);
366  for(int ism = 0; ism < reco.fEMCALL1PhaseInTimeRecalibration->GetEntries(); ism++) fEMCALL1PhaseInTimeRecalibration->AddAt(reco.fEMCALL1PhaseInTimeRecalibration->At(ism), ism);
367  }
368 
369  return *this;
370 }
371 
375 //_____________________________________
377 {
379  {
381  }
382 
384  {
386  }
387 
389  {
391  }
392 
393  if (fEMCALBadChannelMap)
394  {
395  delete fEMCALBadChannelMap;
396  }
397 
398  delete fMatchedTrackIndex ;
399  delete fMatchedClusterIndex ;
400  delete fResidualEta ;
401  delete fResidualPhi ;
402  delete fPIDUtils ;
403 
404  if (fGlobalEventID){
405  delete fGlobalEventID;
406  }
407 
408  InitTrackCuts();
409 }
410 
423 //_______________________________________________________________________________
425  Float_t & amp, Double_t & time,
426  AliVCaloCells* cells)
427 {
428  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
429 
430  if(!geom)
431  {
432  AliError("No instance of the geometry is available");
433  return kFALSE;
434  }
435 
436  if ( absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules() )
437  return kFALSE;
438 
439  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1, status=0;
440 
441  if (!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
442  {
443  // cell absID does not exist
444  amp=0; time = 1.e9;
445  return kFALSE;
446  }
447 
448  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
449 
450  // Do not include bad channels found in analysis,
452  {
453  Bool_t bad = kFALSE;
454 
455  if(fUse1Dmap)
456  bad = GetEMCALChannelStatus1D(absID,status);
457  else
458  bad = GetEMCALChannelStatus(imod, ieta, iphi,status);
459 
460  if ( status > 0 )
461  AliDebug(1,Form("Channel absId %d, status %d, set as bad %d",absID, status, bad));
462 
463  if ( bad ) return kFALSE;
464  }
465  Bool_t isLowGain = !(cells->GetCellHighGain(absID));//HG = false -> LG = true
466 
467  //Recalibrate energy
468  amp = cells->GetCellAmplitude(absID);
470  if(fUse1Drecalib)
472  else
473  amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
474 
475  if(fUseShaperNonlin && isLowGain){
477  }
478  }
479  // Recalibrate time
480  time = cells->GetCellTime(absID);
481  time-=fConstantTimeShift*1e-9; // only in case of old Run1 simulation
482 
483  RecalibrateCellTime(absID,bc,time,isLowGain);
484 
485  //Recalibrate time with L1 phase
487 
488  return kTRUE;
489 }
490 
500 //_____________________________________________________________________________
502  const AliVCluster* cluster,
503  AliVCaloCells* cells)
504 {
505  if (!cluster)
506  {
507  AliInfo("Cluster pointer null!");
508  return kFALSE;
509  }
510 
511  // If the distance to the border is 0 or negative just exit accept all clusters
512  if (cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 )
513  return kTRUE;
514 
515  Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
516  Bool_t shared = kFALSE;
517  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
518 
519  AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
520  absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
521 
522  if (absIdMax==-1) return kFALSE;
523 
524  // Check if the cell is close to the borders:
525  Bool_t okrow = kFALSE;
526  Bool_t okcol = kFALSE;
527 
528  if (iSM < 0 || iphi < 0 || ieta < 0 )
529  {
530  AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
531  iSM,ieta,iphi));
532  return kFALSE; // trick coverity
533  }
534 
535  // Check rows/phi
536  Int_t iPhiLast = 24;
537  if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_Half ) iPhiLast /= 2;
538  else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_3rd ) iPhiLast /= 3;// 1/3 sm case
539  else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Ext ) iPhiLast /= 3;// 1/3 sm case
540 
541  if(iphi >= fNCellsFromEMCALBorder && iphi < iPhiLast - fNCellsFromEMCALBorder) okrow = kTRUE;
542 
543  // Check columns/eta
544  Int_t iEtaLast = 48;
545  if ( !fNoEMCALBorderAtEta0 || geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Standard )
546  {
547  // consider inner border
548  if ( geom->IsDCALSM(iSM) ) iEtaLast = iEtaLast*2/3;
549 
550  if ( ieta > fNCellsFromEMCALBorder && ieta < iEtaLast-fNCellsFromEMCALBorder ) okcol = kTRUE;
551  }
552  else
553  {
554  if (iSM%2==0)
555  {
556  if (ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
557  }
558  else
559  {
560  if (ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE;
561  }
562  }//eta 0 not checked
563 
564  AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
565  fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
566 
567  if (okcol && okrow)
568  {
569  return kTRUE;
570  }
571  else
572  {
573  AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
574 
575  return kFALSE;
576  }
577 }
578 
589 //_______________________________________________________________________________
591  const UShort_t* cellList,
592  Int_t nCells)
593 {
594  if (!fRemoveBadChannels) return kFALSE;
595  if (!fEMCALBadChannelMap) return kFALSE;
596 
597  Int_t icol = -1;
598  Int_t irow = -1;
599  Int_t imod = -1;
600  for (Int_t iCell = 0; iCell<nCells; iCell++)
601  {
602  //Get the column and row
603  Int_t iTower = -1, iIphi = -1, iIeta = -1;
604  geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
605 
606  if (fEMCALBadChannelMap->GetEntries() <= imod) continue;
607 
608  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
609 
610  Int_t status = 0;
611 
612  if(fUse1Dmap){
613  if (GetEMCALChannelStatus1D(cellList[iCell], status))
614  {
615  AliDebug(2,Form("Cluster with bad channel: ID %d, status %d\n",cellList[iCell], status));
616  return kTRUE;
617  }
618  }else{
619  if (GetEMCALChannelStatus(imod, icol, irow, status))
620  {
621  AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d, status %d\n",imod, icol, irow, status));
622  return kTRUE;
623  }
624  }
625  }// cell cluster loop
626 
627  return kFALSE;
628 }
629 
644 //___________________________________________________________________________
646  AliVCaloCells* cells, Int_t bc,
647  Float_t cellMinEn, Bool_t useWeight, Float_t energy )
648 {
649  AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
650 
651  if(!geom)
652  {
653  AliError("No instance of the geometry is available");
654  return -1;
655  }
656 
657  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
658  geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
659  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
660 
661  // Get close cells index, energy and time, not in corners
662 
663  Int_t absID1 = -1;
664  Int_t absID2 = -1;
665 
666  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
667  if ( iphi > 0 ) absID2 = geom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
668 
669  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
670 
671  Int_t absID3 = -1;
672  Int_t absID4 = -1;
673 
674  if ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) )
675  {
676  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
677  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
678  }
679  else if ( ieta == 0 && imod%2 )
680  {
681  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
682  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
683  }
684  else
685  {
686  if ( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
687  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
688  if ( ieta > 0 )
689  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
690  }
691 
692  //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4);
693 
694  Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
695  Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
696 
697  AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
698  AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
699  AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
700  AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
701 
702  if (TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
703  if (TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
704  if (TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
705  if (TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
706 
707  Float_t w1 = 1, w2 = 1, w3 = 1, w4 = 1;
708  if ( useWeight )
709  {
710  w1 = GetCellWeight(ecell1,energy);
711  w2 = GetCellWeight(ecell2,energy);
712  w3 = GetCellWeight(ecell3,energy);
713  w4 = GetCellWeight(ecell4,energy);
714  }
715 
716  if ( ecell1 < cellMinEn || w1 <= 0 ) ecell1 = 0 ;
717  if ( ecell2 < cellMinEn || w2 <= 0 ) ecell2 = 0 ;
718  if ( ecell3 < cellMinEn || w3 <= 0 ) ecell3 = 0 ;
719  if ( ecell4 < cellMinEn || w4 <= 0 ) ecell4 = 0 ;
720 
721  return ecell1+ecell2+ecell3+ecell4;
722 }
723 
724 //________________________________________________________________________________________
734 //________________________________________________________________________________________
736  Int_t & rowDiff, Int_t & colDiff) const
737 {
738  AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
739 
740  if ( !geom )
741  {
742  AliError("No instance of the geometry is available");
743  return -1;
744  }
745 
746  rowDiff = -100;
747  colDiff = -100;
748 
749  if(absId1 == absId2) return kFALSE;
750 
751  // Check if in same SM, if not for sure not same TCard
752  Int_t sm1 = geom->GetSuperModuleNumber(absId1);
753  Int_t sm2 = geom->GetSuperModuleNumber(absId2);
754  if ( sm1 != sm2 ) return kFALSE ;
755 
756  // Get the column and row of each absId
757  Int_t iTower = -1, iIphi = -1, iIeta = -1;
758 
759  Int_t col1, row1;
760  geom->GetCellIndex(absId1,sm1,iTower,iIphi,iIeta);
761  geom->GetCellPhiEtaIndexInSModule(sm1,iTower,iIphi, iIeta,row1,col1);
762 
763  Int_t col2, row2;
764  geom->GetCellIndex(absId2,sm2,iTower,iIphi,iIeta);
765  geom->GetCellPhiEtaIndexInSModule(sm2,iTower,iIphi, iIeta,row2,col2);
766 
767  Int_t row0 = Int_t(row1-row1%8);
768  Int_t col0 = Int_t(col1-col1%2);
769 
770  Int_t rowDiff0 = row2-row0;
771  Int_t colDiff0 = col2-col0;
772 
773  rowDiff = row1-row2;
774  colDiff = col1-col2;
775 
776  // TCard is made by 2x8 towers
777  if ( colDiff0 >=0 && colDiff0 < 2 && rowDiff0 >=0 && rowDiff0 < 8 )
778  {
779 
780  // printf("\t absId (%d,%d), sm %d; col (%d,%d), colDiff %d; row (%d,%d),rowDiff %d\n",
781  // absId1 , absId2, sm1,
782  // col1, col2, colDiff,
783  // row1, row2, rowDiff);
784  return kTRUE ;
785  }
786  else
787  return kFALSE;
788 }
789 
790 //________________________________________________________________________________________
804 //________________________________________________________________________________________
806 (AliVCluster* clus, Int_t absIdMax, AliVCaloCells* cells,
807  Int_t & nDiff, Int_t & nSame,
808  Float_t & eDiff, Float_t & eSame,
809  Float_t emin)
810 {
811  Int_t nCaloCellsPerCluster = clus->GetNCells();
812 
813  nDiff = 0;
814  nSame = 0;
815  Int_t absId = -1;
816  Float_t amp = 0;
817  Int_t rowDiff = -100, colDiff = -100;
818 
819  // Loop on cluster cells count those in same or different T-Card
820  // with respect highest energy cell.
821  for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
822  {
823  absId = clus->GetCellsAbsId()[ipos];
824 
825  amp = cells->GetCellAmplitude(absId);
826 
827  if ( absId == absIdMax || amp < emin ) continue;
828 
829  if ( IsAbsIDsFromTCard(absIdMax,absId,rowDiff,colDiff) )
830  {
831  nSame++;
832  eSame+=amp;
833  }
834  else
835  {
836  nDiff++;
837  eDiff+= amp;
838  }
839 
840  } // cell cluster loop
841 
842 }
843 
854 //_____________________________________________________________________________________________
855 Bool_t AliEMCALRecoUtils::IsExoticCell(Int_t absID, AliVCaloCells* cells, Int_t bc)
856 {
857  if (!fRejectExoticCells) return kFALSE;
858 
859  Float_t ecell = 0;
860  Double_t tcell = 0;
861  Bool_t accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
862 
863  if (!accept) return kTRUE; // reject this cell
864 
865  if (ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
866 
867  Float_t eCross = GetECross(absID,tcell,cells,bc);
868 
869  if (1-eCross/ecell > fExoticCellFraction)
870  {
871  AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",
872  absID,ecell,eCross,1-eCross/ecell));
873  return kTRUE;
874  }
875 
876  return kFALSE;
877 }
878 
888 //___________________________________________________________________
889 Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster,
890  AliVCaloCells *cells,
891  Int_t bc)
892 {
893  if (!cluster)
894  {
895  AliInfo("Cluster pointer null!");
896  return kFALSE;
897  }
898 
899  if (!fRejectExoticCluster) return kFALSE;
900 
901  // Get highest energy tower
902  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
903 
904  if(!geom)
905  {
906  AliError("No instance of the geometry is available");
907  return kFALSE;
908  }
909 
910  Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
911  Bool_t shared = kFALSE;
912  GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
913 
914  return IsExoticCell(absId,cells,bc);
915 }
916 
925 //_______________________________________________________________________
926 Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster)
927 {
928  if (!cluster)
929  {
930  AliInfo("Cluster pointer null!");
931  return 0;
932  }
933 
934  Float_t energy = cluster->E() ;
935  Float_t rdmEnergy = energy ;
936 
937  if (fSmearClusterEnergy)
938  {
939  rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
940  fSmearClusterParam[1] * energy +
941  fSmearClusterParam[2] );
942  AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
943  }
944 
945  return rdmEnergy;
946 }
947 
956 //____________________________________________________________________________
958 {
959  if (!cluster)
960  {
961  AliInfo("Cluster pointer null!");
962  return 0;
963  }
964 
965  Float_t energy = cluster->E();
966  if (energy < 0.100)
967  {
968  // Clusters with less than 100 MeV or negative are not possible
969  AliInfo(Form("Too low cluster energy!, E = %f < 0.100 GeV",energy));
970  return 0;
971  }
972  else if(energy > 300.)
973  {
974  AliInfo(Form("Too high cluster energy!, E = %f GeV, do not apply non linearity",energy));
975  return energy;
976  }
977 
978  switch (fNonLinearityFunction)
979  {
980  case kPi0MC:
981  {
982  //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
983  //fNonLinearityParams[0] = 1.014;
984  //fNonLinearityParams[1] =-0.03329;
985  //fNonLinearityParams[2] =-0.3853;
986  //fNonLinearityParams[3] = 0.5423;
987  //fNonLinearityParams[4] =-0.4335;
988  energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
989  ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
990  exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
991  break;
992  }
993 
994  case kPi0MCv2:
995  {
996  //Non-Linearity correction (from MC with function [0]/((x+[1])^[2]))+1;
997  //fNonLinearityParams[0] = 3.11111e-02;
998  //fNonLinearityParams[1] =-5.71666e-02;
999  //fNonLinearityParams[2] = 5.67995e-01;
1000 
1001  energy *= fNonLinearityParams[0]/TMath::Power(energy+fNonLinearityParams[1],fNonLinearityParams[2])+1;
1002  break;
1003  }
1004 
1005  case kPi0MCv3:
1006  {
1007  //Same as beam test corrected, change parameters
1008  //fNonLinearityParams[0] = 9.81039e-01
1009  //fNonLinearityParams[1] = 1.13508e-01;
1010  //fNonLinearityParams[2] = 1.00173e+00;
1011  //fNonLinearityParams[3] = 9.67998e-02;
1012  //fNonLinearityParams[4] = 2.19381e+02;
1013  //fNonLinearityParams[5] = 6.31604e+01;
1014  //fNonLinearityParams[6] = 1;
1015  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1016 
1017  break;
1018  }
1019 
1020 
1021  case kPi0GammaGamma:
1022  {
1023  //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
1024  //fNonLinearityParams[0] = 1.04;
1025  //fNonLinearityParams[1] = -0.1445;
1026  //fNonLinearityParams[2] = 1.046;
1027  energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
1028  break;
1029  }
1030 
1031  case kPi0GammaConversion:
1032  {
1033  //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
1034  //fNonLinearityParams[0] = 0.139393/0.1349766;
1035  //fNonLinearityParams[1] = 0.0566186;
1036  //fNonLinearityParams[2] = 0.982133;
1037  energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
1038 
1039  break;
1040  }
1041 
1042  case kBeamTest:
1043  {
1044  //From beam test, Alexei's results, for different ZS thresholds
1045  // th=30 MeV; th = 45 MeV; th = 75 MeV
1046  //fNonLinearityParams[0] = 1.007; 1.003; 1.002
1047  //fNonLinearityParams[1] = 0.894; 0.719; 0.797
1048  //fNonLinearityParams[2] = 0.246; 0.334; 0.358
1049  //Rescale the param[0] with 1.03
1050  energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
1051 
1052  break;
1053  }
1054 
1055  case kBeamTestCorrected:
1056  {
1057  // From beam test, corrected for material between beam and EMCAL
1058  //fNonLinearityParams[0] = 0.99078
1059  //fNonLinearityParams[1] = 0.161499;
1060  //fNonLinearityParams[2] = 0.655166;
1061  //fNonLinearityParams[3] = 0.134101;
1062  //fNonLinearityParams[4] = 163.282;
1063  //fNonLinearityParams[5] = 23.6904;
1064  //fNonLinearityParams[6] = 0.978;
1065  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1066 
1067  break;
1068  }
1069 
1070  case kBeamTestCorrectedv2:
1071  {
1072  // From beam test, corrected for material between beam and EMCAL
1073  // Different parametrization to kBeamTestCorrected
1074  //fNonLinearityParams[0] = 0.983504;
1075  //fNonLinearityParams[1] = 0.210106;
1076  //fNonLinearityParams[2] = 0.897274;
1077  //fNonLinearityParams[3] = 0.0829064;
1078  //fNonLinearityParams[4] = 152.299;
1079  //fNonLinearityParams[5] = 31.5028;
1080  //fNonLinearityParams[6] = 0.968;
1081  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1082 
1083  break;
1084  }
1085 
1086  case kBeamTestCorrectedv3:
1087  {
1088  // Same function as kBeamTestCorrected, different default parametrization.
1089  //fNonLinearityParams[0] = 0.976941;
1090  //fNonLinearityParams[1] = 0.162310;
1091  //fNonLinearityParams[2] = 1.08689;
1092  //fNonLinearityParams[3] = 0.0819592;
1093  //fNonLinearityParams[4] = 152.338;
1094  //fNonLinearityParams[5] = 30.9594;
1095  //fNonLinearityParams[6] = 0.9615;
1096  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1097 
1098  break;
1099  }
1100 
1101  case kBeamTestCorrectedv4:
1102  {
1103  // New parametrization of kBeamTestCorrected,
1104  // fitting new points for E>100 GeV.
1105  // I should have same performance as v3 in the low energies
1106  // See EMCal meeting 21/09/2018 slides
1107  // https://indico.cern.ch/event/759154/contributions/3148448/attachments/1721042/2778585/nonLinearityUpdate.pdf
1108  // and jira ticket EMCAL-190
1109  // Not very smart copy pasting the same function for each new parametrization, need to think how to do it better.
1110 
1111 // fNonLinearityParams[0] = 0.9892;
1112 // fNonLinearityParams[1] = 0.1976;
1113 // fNonLinearityParams[2] = 0.865;
1114 // fNonLinearityParams[3] = 0.06775;
1115 // fNonLinearityParams[4] = 156.6;
1116 // fNonLinearityParams[5] = 47.18;
1117 // fNonLinearityParams[6] = 0.97;
1118 
1119  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1120 
1121  break;
1122  }
1123  case kBeamTestNS:
1124  {
1125  // New parametrization of testbeam data points,
1126  // includes also points for E>100 GeV.
1127  // See EMCal meeting 07/12/2018 slides
1128  // https://indico.cern.ch/event/761682/contributions/3245317/attachments/1767706/2870846/2018_12_pp5TeV_NonlinearityStudies_update.pdf
1129 
1130 // fNonLinearityParams[0] = 0.986154;
1131 // fNonLinearityParams[1] = 0.214860;
1132 // fNonLinearityParams[2] = 0.717724;
1133 // fNonLinearityParams[3] = 0.069200;
1134 // fNonLinearityParams[4] = 155.497605;
1135 // fNonLinearityParams[5] = 48.868069;
1136 // fNonLinearityParams[6] = 0.972947;
1137  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1138 
1139  break;
1140  }
1141 
1142  case kSDMv5:
1143  {
1144  //Based on fit to the MC/data using kNoCorrection on the data - utilizes symmetric decay method and kPi0MCv5(MC) - 28 Oct 2013
1145  //fNonLinearityParams[0] = 1.0;
1146  //fNonLinearityParams[1] = 6.64778e-02;
1147  //fNonLinearityParams[2] = 1.570;
1148  //fNonLinearityParams[3] = 9.67998e-02;
1149  //fNonLinearityParams[4] = 2.19381e+02;
1150  //fNonLinearityParams[5] = 6.31604e+01;
1151  //fNonLinearityParams[6] = 1.01286;
1152  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));
1153 
1154  break;
1155  }
1156 
1157  case kPi0MCv5:
1158  {
1159  //Based on comparing MC truth information to the reconstructed energy of clusters.
1160  //fNonLinearityParams[0] = 1.0;
1161  //fNonLinearityParams[1] = 6.64778e-02;
1162  //fNonLinearityParams[2] = 1.570;
1163  //fNonLinearityParams[3] = 9.67998e-02;
1164  //fNonLinearityParams[4] = 2.19381e+02;
1165  //fNonLinearityParams[5] = 6.31604e+01;
1166  //fNonLinearityParams[6] = 1.01286;
1167  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1168 
1169  break;
1170  }
1171 
1172  case kSDMv6:
1173  {
1174  //Based on fit to the MC/data using kNoCorrection on the data
1175  // - utilizes symmetric decay method and kPi0MCv6(MC) - 09 Dec 2014
1176  // - parameters constrained by the test beam data as well
1177  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1178  //fNonLinearityParams[0] = 1.0;
1179  //fNonLinearityParams[1] = 0.237767;
1180  //fNonLinearityParams[2] = 0.651203;
1181  //fNonLinearityParams[3] = 0.183741;
1182  //fNonLinearityParams[4] = 155.427;
1183  //fNonLinearityParams[5] = 17.0335;
1184  //fNonLinearityParams[6] = 0.987054;
1185  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1186 
1187  break;
1188  }
1189 
1190  case kPi0MCv6:
1191  {
1192  //Based on comparing MC truth information to the reconstructed energy of clusters.
1193  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1194  //fNonLinearityParams[0] = 1.0;
1195  //fNonLinearityParams[1] = 0.0797873;
1196  //fNonLinearityParams[2] = 1.68322;
1197  //fNonLinearityParams[3] = 0.0806098;
1198  //fNonLinearityParams[4] = 244.586;
1199  //fNonLinearityParams[5] = 116.938;
1200  //fNonLinearityParams[6] = 1.00437;
1201  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1202 
1203  break;
1204  }
1205 
1206  case kPi0MCNS:
1207  {
1208  // New parametrization of testbeam MC points,
1209  // includes also points for E>100 GeV.
1210  // See EMCal meeting 07/12/2018 slides
1211  // https://indico.cern.ch/event/761682/contributions/3245317/attachments/1767706/2870846/2018_12_pp5TeV_NonlinearityStudies_update.pdf
1212  //fNonLinearityParams[0] = 1.009121;
1213  //fNonLinearityParams[1] = 0.083153;
1214  //fNonLinearityParams[2] = 1.444362;
1215  //fNonLinearityParams[3] = 0.100294;
1216  //fNonLinearityParams[4] = 416.897753;
1217  //fNonLinearityParams[5] = 324.246101;
1218  //fNonLinearityParams[6] = 1.004055;
1219  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1220 
1221  break;
1222  }
1223 
1224  case kPCMv1:
1225  {
1226  //based on symmetric decays of pi0 meson
1227  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1228  // parameters vary from MC to MC
1229  //fNonLinearityParams[0] = 0.984876;
1230  //fNonLinearityParams[1] = -9.999609;
1231  //fNonLinearityParams[2] = -4.999891;
1232  //fNonLinearityParams[3] = 0.;
1233  //fNonLinearityParams[4] = 0.;
1234  //fNonLinearityParams[5] = 0.;
1235  //fNonLinearityParams[6] = 0.;
1236  energy /= TMath::Power(fNonLinearityParams[0] + exp(fNonLinearityParams[1] + fNonLinearityParams[2]*energy),2);//result coming from calo-conv method
1237 
1238  break;
1239  }
1240 
1241  case kPCMplusBTCv1:
1242  {
1243  //convolution of TestBeamCorrectedv3 with PCM method
1244  //Based on comparing MC truth information to the reconstructed energy of clusters.
1245  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1246  // parameters vary from MC to MC
1247  //fNonLinearityParams[0] = 0.976941;
1248  //fNonLinearityParams[1] = 0.162310;
1249  //fNonLinearityParams[2] = 1.08689;
1250  //fNonLinearityParams[3] = 0.0819592;
1251  //fNonLinearityParams[4] = 152.338;
1252  //fNonLinearityParams[5] = 30.9594;
1253  //fNonLinearityParams[6] = 0.9615;
1254  //fNonLinearityParams[7] = 0.984876;
1255  //fNonLinearityParams[8] = -9.999609;
1256  //fNonLinearityParams[9] = -4.999891;
1257  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1258  energy /= TMath::Power(fNonLinearityParams[7] + exp(fNonLinearityParams[8] + fNonLinearityParams[9]*energy),2);//result coming from calo-conv method
1259 
1260  break;
1261  }
1262 
1263  case kPCMsysv1:
1264  {
1265  // Systematic variation of kPCMv1
1266  //Based on comparing MC truth information to the reconstructed energy of clusters.
1267  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1268  // parameters vary from MC to MC
1269  //fNonLinearityParams[0] = 0.0;
1270  //fNonLinearityParams[1] = 1.0;
1271  //fNonLinearityParams[2] = 1.0;
1272  //fNonLinearityParams[3] = 0.0;
1273  //fNonLinearityParams[4] = 1.0;
1274  //fNonLinearityParams[5] = 0.0;
1275  //fNonLinearityParams[6] = 0.0;
1276  energy /= TMath::Power( (fNonLinearityParams[0] + fNonLinearityParams[1] * TMath::Power(energy,fNonLinearityParams[2]) ) /
1277  (fNonLinearityParams[3] + fNonLinearityParams[4] * TMath::Power(energy,fNonLinearityParams[5]) ) + fNonLinearityParams[6], 2);//result coming from calo-conv method
1278 
1279  break;
1280  }
1281 
1282 
1283  case kNoCorrection:
1284  AliDebug(2,"No correction on the energy\n");
1285  break;
1286 
1287  }
1288 
1289  return energy;
1290 }
1291 
1296 //__________________________________________________
1298 {
1299  if (fNonLinearityFunction == kPi0MC) {
1300  fNonLinearityParams[0] = 1.014;
1301  fNonLinearityParams[1] = -0.03329;
1302  fNonLinearityParams[2] = -0.3853;
1303  fNonLinearityParams[3] = 0.5423;
1304  fNonLinearityParams[4] = -0.4335;
1305  }
1306 
1308  fNonLinearityParams[0] = 3.11111e-02;
1309  fNonLinearityParams[1] =-5.71666e-02;
1310  fNonLinearityParams[2] = 5.67995e-01;
1311  }
1312 
1314  fNonLinearityParams[0] = 9.81039e-01;
1315  fNonLinearityParams[1] = 1.13508e-01;
1316  fNonLinearityParams[2] = 1.00173e+00;
1317  fNonLinearityParams[3] = 9.67998e-02;
1318  fNonLinearityParams[4] = 2.19381e+02;
1319  fNonLinearityParams[5] = 6.31604e+01;
1320  fNonLinearityParams[6] = 1;
1321  }
1322 
1324  fNonLinearityParams[0] = 1.04;
1325  fNonLinearityParams[1] = -0.1445;
1326  fNonLinearityParams[2] = 1.046;
1327  }
1328 
1330  fNonLinearityParams[0] = 0.139393;
1331  fNonLinearityParams[1] = 0.0566186;
1332  fNonLinearityParams[2] = 0.982133;
1333  }
1334 
1336  if (fNonLinearThreshold == 30) {
1337  fNonLinearityParams[0] = 1.007;
1338  fNonLinearityParams[1] = 0.894;
1339  fNonLinearityParams[2] = 0.246;
1340  }
1341  if (fNonLinearThreshold == 45) {
1342  fNonLinearityParams[0] = 1.003;
1343  fNonLinearityParams[1] = 0.719;
1344  fNonLinearityParams[2] = 0.334;
1345  }
1346  if (fNonLinearThreshold == 75) {
1347  fNonLinearityParams[0] = 1.002;
1348  fNonLinearityParams[1] = 0.797;
1349  fNonLinearityParams[2] = 0.358;
1350  }
1351  }
1352 
1354  fNonLinearityParams[0] = 0.99078;
1355  fNonLinearityParams[1] = 0.161499;
1356  fNonLinearityParams[2] = 0.655166;
1357  fNonLinearityParams[3] = 0.134101;
1358  fNonLinearityParams[4] = 163.282;
1359  fNonLinearityParams[5] = 23.6904;
1360  fNonLinearityParams[6] = 0.978;
1361  }
1362 
1364  // Parameters until November 2015, use now kBeamTestCorrectedv3
1365  fNonLinearityParams[0] = 0.983504;
1366  fNonLinearityParams[1] = 0.210106;
1367  fNonLinearityParams[2] = 0.897274;
1368  fNonLinearityParams[3] = 0.0829064;
1369  fNonLinearityParams[4] = 152.299;
1370  fNonLinearityParams[5] = 31.5028;
1371  fNonLinearityParams[6] = 0.968;
1372  }
1373 
1375 
1376  // New parametrization of kBeamTestCorrected
1377  // excluding point at 0.5 GeV from Beam Test Data
1378  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1379 
1380  fNonLinearityParams[0] = 0.976941;
1381  fNonLinearityParams[1] = 0.162310;
1382  fNonLinearityParams[2] = 1.08689;
1383  fNonLinearityParams[3] = 0.0819592;
1384  fNonLinearityParams[4] = 152.338;
1385  fNonLinearityParams[5] = 30.9594;
1386  fNonLinearityParams[6] = 0.9615;
1387  }
1388 
1390 
1391  // New parametrization of kBeamTestCorrected,
1392  // fitting new points for E>100 GeV.
1393  // I should have same performance as v3 in the low energies
1394  // See EMCal meeting 21/09/2018 slides
1395  // https://indico.cern.ch/event/759154/contributions/3148448/attachments/1721042/2778585/nonLinearityUpdate.pdf
1396  // and jira ticket EMCAL-190
1397 
1398  fNonLinearityParams[0] = 0.9892;
1399  fNonLinearityParams[1] = 0.1976;
1400  fNonLinearityParams[2] = 0.865;
1401  fNonLinearityParams[3] = 0.06775;
1402  fNonLinearityParams[4] = 156.6;
1403  fNonLinearityParams[5] = 47.18;
1404  fNonLinearityParams[6] = 0.97;
1405  }
1406 
1408 
1409  // New parametrization of testbeam data points,
1410  // includes also points for E>100 GeV.
1411  // See EMCal meeting 07/12/2018 slides
1412  // https://indico.cern.ch/event/761682/contributions/3245317/attachments/1767706/2870846/2018_12_pp5TeV_NonlinearityStudies_update.pdf
1413 
1414  fNonLinearityParams[0] = 0.986154;
1415  fNonLinearityParams[1] = 0.214860;
1416  fNonLinearityParams[2] = 0.717724;
1417  fNonLinearityParams[3] = 0.069200;
1418  fNonLinearityParams[4] = 155.497605;
1419  fNonLinearityParams[5] = 48.868069;
1420  fNonLinearityParams[6] = 0.972947;
1421  }
1422 
1423  if (fNonLinearityFunction == kSDMv5) {
1424  fNonLinearityParams[0] = 1.0;
1425  fNonLinearityParams[1] = 6.64778e-02;
1426  fNonLinearityParams[2] = 1.570;
1427  fNonLinearityParams[3] = 9.67998e-02;
1428  fNonLinearityParams[4] = 2.19381e+02;
1429  fNonLinearityParams[5] = 6.31604e+01;
1430  fNonLinearityParams[6] = 1.01286;
1431  }
1432 
1434  fNonLinearityParams[0] = 1.0;
1435  fNonLinearityParams[1] = 6.64778e-02;
1436  fNonLinearityParams[2] = 1.570;
1437  fNonLinearityParams[3] = 9.67998e-02;
1438  fNonLinearityParams[4] = 2.19381e+02;
1439  fNonLinearityParams[5] = 6.31604e+01;
1440  fNonLinearityParams[6] = 1.01286;
1441  }
1442 
1443  if (fNonLinearityFunction == kSDMv6) {
1444  fNonLinearityParams[0] = 1.0;
1445  fNonLinearityParams[1] = 0.237767;
1446  fNonLinearityParams[2] = 0.651203;
1447  fNonLinearityParams[3] = 0.183741;
1448  fNonLinearityParams[4] = 155.427;
1449  fNonLinearityParams[5] = 17.0335;
1450  fNonLinearityParams[6] = 0.987054;
1451  }
1452 
1454  fNonLinearityParams[0] = 1.0;
1455  fNonLinearityParams[1] = 0.0797873;
1456  fNonLinearityParams[2] = 1.68322;
1457  fNonLinearityParams[3] = 0.0806098;
1458  fNonLinearityParams[4] = 244.586;
1459  fNonLinearityParams[5] = 116.938;
1460  fNonLinearityParams[6] = 1.00437;
1461  }
1462 
1464  fNonLinearityParams[0] = 1.009121;
1465  fNonLinearityParams[1] = 0.083153;
1466  fNonLinearityParams[2] = 1.444362;
1467  fNonLinearityParams[3] = 0.100294;
1468  fNonLinearityParams[4] = 416.897753;
1469  fNonLinearityParams[5] = 324.246101;
1470  fNonLinearityParams[6] = 1.004055;
1471  }
1472 
1473 if (fNonLinearityFunction == kPCMv1) {
1474  //parameters change from MC production to MC production, they need to set for each period
1475  fNonLinearityParams[0] = 0.984876;
1476  fNonLinearityParams[1] = -9.999609;
1477  fNonLinearityParams[2] = -4.999891;
1478  fNonLinearityParams[3] = 0.;
1479  fNonLinearityParams[4] = 0.;
1480  fNonLinearityParams[5] = 0.;
1481  fNonLinearityParams[6] = 0.;
1482  }
1483 
1485  // test beam corrected values convoluted with symmetric meson decays values
1486  // for test beam:
1487  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1488  // for PCM method:
1489  // https://aliceinfo.cern.ch/Notes/node/211
1490  fNonLinearityParams[0] = 0.976941;
1491  fNonLinearityParams[1] = 0.162310;
1492  fNonLinearityParams[2] = 1.08689;
1493  fNonLinearityParams[3] = 0.0819592;
1494  fNonLinearityParams[4] = 152.338;
1495  fNonLinearityParams[5] = 30.9594;
1496  fNonLinearityParams[6] = 0.9615;
1497  fNonLinearityParams[7] = 0.984876;
1498  fNonLinearityParams[8] = -9.999609;
1499  fNonLinearityParams[9] = -4.999891;
1500  }
1501 
1503  //systematics for kPCMv1
1504  // for PCM method:
1505  // https://aliceinfo.cern.ch/Notes/node/211
1506  fNonLinearityParams[0] = 0.0;
1507  fNonLinearityParams[1] = 1.0;
1508  fNonLinearityParams[2] = 1.0;
1509  fNonLinearityParams[3] = 0.0;
1510  fNonLinearityParams[4] = 1.0;
1511  fNonLinearityParams[5] = 0.0;
1512  fNonLinearityParams[6] = 0.0;
1513  }
1514 }
1515 
1523 //____________________________________________________________________
1525 {
1526  fBadStatusSelection[0] = all; // declare all as bad if true, never mind the other settings
1527  fBadStatusSelection[1] = dead;
1528  fBadStatusSelection[2] = hot;
1529  fBadStatusSelection[3] = warm;
1530 }
1531 
1542 //____________________________________________________________________
1544 {
1545  if(!fUse1Dmap){
1546  if(fEMCALBadChannelMap)
1547  status = (Int_t) ((TH2I*)fEMCALBadChannelMap->At(iSM))->GetBinContent(iCol,iRow);
1548  else
1549  status = 0; // Channel is ok by default
1550  }else{
1551  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
1552  Int_t CellID = geom->GetAbsCellIdFromCellIndexes(iSM, iRow, iCol);
1553  if(fEMCALBadChannelMap)
1554  status = (Int_t) ((TH1C*)fEMCALBadChannelMap->At(0))->GetBinContent(CellID);
1555  else
1556  status = 0; // Channel is ok by default
1557  }
1558 
1559  if ( status == AliCaloCalibPedestal::kAlive )
1560  {
1561  return kFALSE; // Good channel
1562  }
1563  else
1564  {
1565  if ( fBadStatusSelection[0] == kTRUE )
1566  {
1567  return kTRUE; // consider bad hot, dead and warm
1568  }
1569  else
1570  {
1571  if ( fBadStatusSelection[AliCaloCalibPedestal::kDead] == kTRUE &&
1572  status == AliCaloCalibPedestal::kDead )
1573  return kTRUE; // consider bad dead
1574  else if ( fBadStatusSelection[AliCaloCalibPedestal::kHot] == kTRUE &&
1575  status == AliCaloCalibPedestal::kHot )
1576  return kTRUE; // consider bad hot
1577  else if ( fBadStatusSelection[AliCaloCalibPedestal::kWarning] == kTRUE &&
1578  status == AliCaloCalibPedestal::kWarning )
1579  return kTRUE; // consider bad warm
1580  }
1581  }
1582 
1583  AliWarning(Form("Careful, bad channel selection not properly done: ism %d, icol %d, irow %d, status %d,\n"
1584  " fBadAll %d, fBadHot %d, fBadWarm %d, fBadDead %d",
1585  iSM, iCol, iRow, status,
1588 
1589  return kFALSE; // if everything fails, accept it.
1590 }
1591 
1600 //____________________________________________________________________
1602 {
1603  if(fEMCALBadChannelMap)
1604  status = (Int_t) ((TH1C*)fEMCALBadChannelMap->At(0))->GetBinContent(iCell);
1605  else
1606  status = 0; // Channel is ok by default
1607 
1608  if ( status == AliCaloCalibPedestal::kAlive )
1609  {
1610  return kFALSE; // Good channel
1611  }
1612  else
1613  {
1614  if ( fBadStatusSelection[0] == kTRUE )
1615  {
1616  return kTRUE; // consider bad hot, dead and warm
1617  }
1618  else
1619  {
1620  if ( fBadStatusSelection[AliCaloCalibPedestal::kDead] == kTRUE &&
1621  status == AliCaloCalibPedestal::kDead )
1622  return kTRUE; // consider bad dead
1623  else if ( fBadStatusSelection[AliCaloCalibPedestal::kHot] == kTRUE &&
1624  status == AliCaloCalibPedestal::kHot )
1625  return kTRUE; // consider bad hot
1626  else if ( fBadStatusSelection[AliCaloCalibPedestal::kWarning] == kTRUE &&
1627  status == AliCaloCalibPedestal::kWarning )
1628  return kTRUE; // consider bad warm
1629  }
1630  }
1631 
1632  AliWarning(Form("Careful, bad channel selection not properly done: icell %d, status %d,\n"
1633  " fBadAll %d, fBadHot %d, fBadWarm %d, fBadDead %d",
1634  iCell, status,
1637 
1638  return kFALSE; // if everything fails, accept it.
1639 }
1640 
1654 //____________________________________________________________________
1655 void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
1656  AliVCaloCells* cells,
1657  const AliVCluster* clu,
1658  Int_t & absId,
1659  Int_t & iSupMod,
1660  Int_t & ieta,
1661  Int_t & iphi,
1662  Bool_t & shared)
1663 {
1664  Double_t eMax = -1.;
1665  Double_t eCell = -1.;
1666  Float_t fraction = 1.;
1667  Float_t recalFactor = 1.;
1668  Int_t cellAbsId = -1 ;
1669 
1670  Int_t iTower = -1;
1671  Int_t iIphi = -1;
1672  Int_t iIeta = -1;
1673  Int_t iSupMod0= -1;
1674 
1675  if (!clu)
1676  {
1677  AliInfo("Cluster pointer null!");
1678  absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
1679  return;
1680  }
1681 
1682  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1683  {
1684  cellAbsId = clu->GetCellAbsId(iDig);
1685  fraction = clu->GetCellAmplitudeFraction(iDig);
1686  //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
1687  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1688 
1689  geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1690  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1691 
1692  if (iDig==0)
1693  {
1694  iSupMod0=iSupMod;
1695  } else if (iSupMod0!=iSupMod)
1696  {
1697  shared = kTRUE;
1698  //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
1699  }
1700 
1702  if(fUse1Drecalib)
1703  recalFactor = GetEMCALChannelRecalibrationFactor1D(cellAbsId);
1704  else
1705  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1706  }
1707 
1708  eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
1709  //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
1710  if (eCell > eMax)
1711  {
1712  eMax = eCell;
1713  absId = cellAbsId;
1714  //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
1715  }
1716  }// cell loop
1717 
1718  //Get from the absid the supermodule, tower and eta/phi numbers
1719  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1720  //Gives SuperModule and Tower numbers
1721  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1722  //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
1723  //printf("Max end---\n");
1724 }
1725 
1734 //_________________________________________________________
1736 {
1737  if (eCell > 0 && eCluster > 0)
1738  {
1739  if ( fW0 > 0 ) return TMath::Max( 0., fW0 + TMath::Log( eCell / eCluster ) ) ;
1740  else return TMath::Log( eCluster / eCell ) ;
1741  }
1742  else
1743  return 0. ;
1744 }
1745 
1749 //______________________________________
1751 {
1752  fParticleType = kPhoton;
1753  fPosAlgo = kUnchanged;
1754  fW0 = 4.5;
1755 
1757  fNonLinearThreshold = 30;
1758 
1759  fExoticCellFraction = 0.97;
1760  fExoticCellDiffTime = 1e6;
1762 
1763  fAODFilterMask = 128;
1764  fAODHybridTracks = kFALSE;
1765  fAODTPCOnlyTracks = kTRUE;
1766 
1767  fCutEtaPhiSum = kTRUE;
1768  fCutEtaPhiSeparate = kFALSE;
1769 
1770  fCutR = 0.05;
1771  fCutEta = 0.025;
1772  fCutPhi = 0.05;
1773 
1774  fClusterWindow = 100;
1775  fMass = 0.139;
1776 
1777  fStepSurface = 20.;
1778  fStepCluster = 5.;
1780 
1781  fCutMinTrackPt = 0;
1782  fCutMinNClusterTPC = -1;
1783  fCutMinNClusterITS = -1;
1784 
1785  fCutMaxChi2PerClusterTPC = 1e10;
1786  fCutMaxChi2PerClusterITS = 1e10;
1787 
1788  fCutRequireTPCRefit = kFALSE;
1789  fCutRequireITSRefit = kFALSE;
1790  fCutAcceptKinkDaughters = kFALSE;
1791 
1792  fCutMaxDCAToVertexXY = 1e10;
1793  fCutMaxDCAToVertexZ = 1e10;
1794  fCutDCAToVertex2D = kFALSE;
1795 
1796  fCutRequireITSStandAlone = kFALSE; //MARCEL
1797  fCutRequireITSpureSA = kFALSE; //Marcel
1798 
1799  // Misalignment matrices
1800  for (Int_t i = 0; i < 15 ; i++)
1801  {
1802  fMisalTransShift[i] = 0.;
1803  fMisalRotShift[i] = 0.;
1804  }
1805 
1806  // Non linearity
1807  for (Int_t i = 0; i < 10 ; i++) fNonLinearityParams[i] = 0.;
1808 
1809  // For kBeamTestCorrectedv2 case, but default is no correction
1810  fNonLinearityParams[0] = 0.983504;
1811  fNonLinearityParams[1] = 0.210106;
1812  fNonLinearityParams[2] = 0.897274;
1813  fNonLinearityParams[3] = 0.0829064;
1814  fNonLinearityParams[4] = 152.299;
1815  fNonLinearityParams[5] = 31.5028;
1816  fNonLinearityParams[6] = 0.968;
1817 
1818  // Cluster energy smearing
1819  fSmearClusterEnergy = kFALSE;
1820  fSmearClusterParam[0] = 0.07; // * sqrt E term
1821  fSmearClusterParam[1] = 0.00; // * E term
1822  fSmearClusterParam[2] = 0.00; // constant
1823 }
1824 
1828 //_____________________________________________________
1830 {
1831  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1832 
1833  // In order to avoid rewriting the same histograms
1834  Bool_t oldStatus = TH1::AddDirectoryStatus();
1835  TH1::AddDirectory(kFALSE);
1836 
1838  for (int i = 0; i < 22; i++)
1839  fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1840  Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
1841  //Init the histograms with 1
1842  for (Int_t sm = 0; sm < 22; sm++)
1843  {
1844  for (Int_t i = 0; i < 48; i++)
1845  {
1846  for (Int_t j = 0; j < 24; j++)
1847  {
1849  }
1850  }
1851  }
1852 
1853  fEMCALRecalibrationFactors->SetOwner(kTRUE);
1854  fEMCALRecalibrationFactors->Compress();
1855 
1856  // In order to avoid rewriting the same histograms
1857  TH1::AddDirectory(oldStatus);
1858 }
1862 //_____________________________________________________
1864 {
1865  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors1D()");
1866 
1867  // In order to avoid rewriting the same histograms
1868  Bool_t oldStatus = TH1::AddDirectoryStatus();
1869  TH1::AddDirectory(kFALSE);
1870 
1872 
1873  fEMCALRecalibrationFactors->Add(new TH1S("EMCALRecalFactors","EMCALRecalFactors", 48*24*22,0.,48*24*22));
1874 
1875  //Init the histograms with 1
1876  for (UInt_t icell = 0; icell < 48*24*22; icell++)
1877  {
1879  }
1880 
1881  fEMCALRecalibrationFactors->SetOwner(kTRUE);
1882  fEMCALRecalibrationFactors->Compress();
1883 
1884  // In order to avoid rewriting the same histograms
1885  TH1::AddDirectory(oldStatus);
1886 }
1887 
1891 //_________________________________________________________
1893 {
1894  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1895 
1896  // In order to avoid rewriting the same histograms
1897  Bool_t oldStatus = TH1::AddDirectoryStatus();
1898  TH1::AddDirectory(kFALSE);
1899 
1900  if(fDoUseMergedBC){
1901 
1904 
1905  fEMCALTimeRecalibrationFactors->Add(new TH1S("hAllTimeAv",
1906  "hAllTimeAv",
1907  48*24*22,0.,48*24*22) );
1908  // Init the histograms with 0
1909  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1910  SetEMCALChannelTimeRecalibrationFactor(0,iCh,0.,kFALSE);
1911 
1912  if(fLowGain) {
1913  fEMCALTimeRecalibrationFactors->Add(new TH1F("hAllTimeAvLG",
1914  "hAllTimeAvLG",
1915  48*24*22,0.,48*24*22) );
1916  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1918  }
1919 
1920  }else{
1923 
1924  for (int i = 0; i < 4; i++)
1925  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1926  Form("hAllTimeAvBC%d",i),
1927  48*24*22,0.,48*24*22) );
1928  // Init the histograms with 0
1929  for (Int_t iBC = 0; iBC < 4; iBC++)
1930  {
1931  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1932  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kFALSE);
1933  }
1934 
1935  if(fLowGain) {
1936  for (int iBC = 0; iBC < 4; iBC++) {
1937  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvLGBC%d",iBC),
1938  Form("hAllTimeAvLGBC%d",iBC),
1939  48*24*22,0.,48*24*22) );
1940  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1941  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kTRUE);
1942  }
1943  }
1944 
1945  }
1946 
1947  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1948  fEMCALTimeRecalibrationFactors->Compress();
1949 
1950  // In order to avoid rewriting the same histograms
1951  TH1::AddDirectory(oldStatus);
1952 }
1953 
1957 //____________________________________________________
1959 {
1960  AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1961 
1962  // In order to avoid rewriting the same histograms
1963  Bool_t oldStatus = TH1::AddDirectoryStatus();
1964  TH1::AddDirectory(kFALSE);
1965 
1966  fEMCALBadChannelMap = new TObjArray(22);
1967  //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
1968 
1969  for (int i = 0; i < 22; i++)
1970  fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1971 
1972  fEMCALBadChannelMap->SetOwner(kTRUE);
1973  fEMCALBadChannelMap->Compress();
1974 
1975  // In order to avoid rewriting the same histograms
1976  TH1::AddDirectory(oldStatus);
1977 }
1978 
1982 //____________________________________________________
1984 {
1985  AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap1D()");
1986 
1987  fUse1Dmap = kTRUE;
1988  // In order to avoid rewriting the same histograms
1989  Bool_t oldStatus = TH1::AddDirectoryStatus();
1990  TH1::AddDirectory(kFALSE);
1991 
1992  fEMCALBadChannelMap = new TObjArray(1);
1993 
1994  fEMCALBadChannelMap->Add(new TH1C("EMCALBadChannelMap","EMCALBadChannelMap", 48*24*22,0.,48*24*22));
1995 
1996  fEMCALBadChannelMap->SetOwner(kTRUE);
1997  fEMCALBadChannelMap->Compress();
1998 
1999  // In order to avoid rewriting the same histograms
2000  TH1::AddDirectory(oldStatus);
2001 }
2002 
2006 //___________________________________________________________
2008 {
2009  AliDebug(2,"AliEMCALRecoUtils::InitEMCALL1PhaseInTimeRecalibrationFactors()");
2010 
2011  // In order to avoid rewriting the same histograms
2012  Bool_t oldStatus = TH1::AddDirectoryStatus();
2013  TH1::AddDirectory(kFALSE);
2014 
2016 
2017  fEMCALL1PhaseInTimeRecalibration->Add(new TH1C("h0","EMCALL1phaseForSM", 22, 0, 22));
2018 
2019  for (Int_t i = 0; i < 22; i++) //loop over SMs, default value = 0
2021 
2022  fEMCALL1PhaseInTimeRecalibration->SetOwner(kTRUE);
2024 
2025  // In order to avoid rewriting the same histograms
2026  TH1::AddDirectory(oldStatus);
2027 }
2028 
2038 //____________________________________________________________________________
2039 void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
2040  AliVCluster * cluster,
2041  AliVCaloCells * cells,
2042  Int_t bc)
2043 {
2044  if (!cluster)
2045  {
2046  AliInfo("Cluster pointer null!");
2047  return;
2048  }
2049 
2050  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
2051  UShort_t * index = cluster->GetCellsAbsId() ;
2052  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
2053  Int_t ncells = cluster->GetNCells();
2054 
2055  // Initialize some used variables
2056  Float_t energy = 0;
2057  Int_t absId =-1;
2058  Int_t icol =-1, irow =-1, imod=1;
2059  Float_t factor = 1, frac = 0;
2060  Int_t absIdMax = -1;
2061  Float_t emax = 0;
2062 
2063  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
2064  for (Int_t icell = 0; icell < ncells; icell++)
2065  {
2066  absId = index[icell];
2067  frac = fraction[icell];
2068  if (frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
2069 
2071  {
2072  // Energy
2073  Int_t iTower = -1, iIphi = -1, iIeta = -1;
2074  geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
2075  if (fEMCALRecalibrationFactors->GetEntries() <= imod)
2076  continue;
2077  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
2078 
2079  if(fUse1Drecalib)
2080  factor = GetEMCALChannelRecalibrationFactor1D(absId);
2081  else
2082  factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
2083 
2084  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
2085  imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
2086 
2087  }
2088 
2089  energy += cells->GetCellAmplitude(absId)*factor*frac;
2090 
2091  if (emax < cells->GetCellAmplitude(absId)*factor*frac)
2092  {
2093  emax = cells->GetCellAmplitude(absId)*factor*frac;
2094  absIdMax = absId;
2095  }
2096  }
2097 
2098  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
2099 
2100  cluster->SetE(energy);
2101 
2102  // Recalculate time of cluster
2103  Double_t timeorg = cluster->GetTOF();
2104  Bool_t isLowGain = !(cells->GetCellHighGain(absIdMax));//HG = false -> LG = true
2105 
2106  Double_t time = cells->GetCellTime(absIdMax);
2108  RecalibrateCellTime(absIdMax,bc,time,isLowGain);
2109  time-=fConstantTimeShift*1e-9; // only in case of Run1 old simulations
2110 
2111  // Recalibrate time with L1 phase
2114 
2115  cluster->SetTOF(time);
2116 
2117  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
2118 }
2119 
2127 //_______________________________________________________________________
2128 void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, Int_t bc)
2129 {
2131  return;
2132 
2133  if (!cells)
2134  {
2135  AliInfo("Cells pointer null!");
2136  return;
2137  }
2138 
2139  Short_t absId =-1;
2140  Bool_t accept = kFALSE;
2141  Float_t ecell = 0;
2142  Double_t tcell = 0;
2143  Double_t ecellin = 0;
2144  Double_t tcellin = 0;
2145  Int_t mclabel = -1;
2146  Double_t efrac = 0;
2147 
2148  Int_t nEMcell = cells->GetNumberOfCells() ;
2149  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
2150  {
2151  cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
2152 
2153  accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
2154  if (!accept)
2155  {
2156  ecell = 0;
2157  tcell = -1;
2158  }
2159 
2160  // Set new values
2161  cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
2162  }
2163 
2164  fCellsRecalibrated = kTRUE;
2165 }
2166 
2173 //_______________________________________________________________________
2175 {
2176  // The following conversion factor needs to be applied to go from energy to ADC
2177  // AliEMCALCalibData::fADCchannelRef = 0.0162;
2178 
2179  if(Emeas<40){
2180  return Emeas*16.3/16;
2181  }
2182  Float_t par[]={1, 29.8279, 0.607704, 0.00164896, -2.28595e-06, -8.54664e-10, 5.50191e-12, -3.28098e-15};
2183  Float_t x = par[0]*Emeas/EcalibHG/16/0.0162;
2184 
2185  Float_t res=par[1];
2186  res+=par[2]*x;
2187  res+=par[3]*x*x;
2188  res+=par[4]*x*x*x;
2189  res+=par[5]*x*x*x*x;
2190  res+=par[6]*x*x*x*x*x;
2191  res+=par[7]*x*x*x*x*x*x;
2192 
2193  return EcalibHG*16.3*res*0.0162;
2194 }
2195 
2204 //_______________________________________________________________________________________________________
2205 void AliEMCALRecoUtils::RecalibrateCellTime(Int_t absId, Int_t bc, Double_t & celltime, Bool_t isLGon) const
2206 {
2207  if (!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0) {
2208  if(fLowGain)
2209  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,isLGon)*1.e-9;
2210  else
2211  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,kFALSE)*1.e-9;
2212  }
2213 }
2214 
2223 //_______________________________________________________________________________________________________
2225 {
2226  if (!fCellsRecalibrated && IsL1PhaseInTimeRecalibrationOn() && bc >= 0)
2227  {
2228  bc=bc%4;
2229 
2230  Float_t offsetPerSM=0.;
2231  Int_t l1PhaseShift = GetEMCALL1PhaseInTimeRecalibrationForSM(iSM,par);
2232  Int_t l1Phase=l1PhaseShift & 3; //bit operation
2233 
2234  if(bc >= l1Phase)
2235  offsetPerSM = (bc - l1Phase)*25;
2236  else
2237  offsetPerSM = (bc - l1Phase + 4)*25;
2238 
2239  Int_t l1shiftOffset=l1PhaseShift>>2; //bit operation
2240  l1shiftOffset*=25;
2241 
2242  celltime -= offsetPerSM*1.e-9;
2243  celltime -= l1shiftOffset*1.e-9;
2244  }
2245 }
2246 
2262 //______________________________________________________________________________
2263 void AliEMCALRecoUtils::RecalculateCellLabelsRemoveAddedGenerator( Int_t absID, AliVCluster* clus, AliMCEvent* mc,
2264  Float_t & amp, TArrayI & labeArr, TArrayF & eDepArr) const
2265 {
2266  TString genName ;
2267  Float_t eDepFrac[4];
2268 
2269  Float_t edepTotFrac = 1;
2270  Bool_t found = kFALSE;
2271  Float_t ampOrg = amp;
2272 
2273  //
2274  // Get the energy deposition fraction from cluster.
2275  //
2276  for(Int_t icluscell = 0; icluscell < clus->GetNCells(); icluscell++ )
2277  {
2278  if(absID == clus->GetCellAbsId(icluscell))
2279  {
2280  clus->GetCellMCEdepFractionArray(icluscell,eDepFrac);
2281 
2282  found = kTRUE;
2283 
2284  break;
2285  }
2286  }
2287 
2288  if ( !found )
2289  {
2290  AliWarning(Form("Cell abs ID %d NOT found in cluster",absID));
2291  return;
2292  }
2293 
2294  //
2295  // Check if there is a particle contribution from a given generator name.
2296  // If it is not one of the selected generators,
2297  // remove the constribution from the cell.
2298  //
2299  if ( mc && fNMCGenerToAccept > 0 )
2300  {
2301  //printf("Accept contribution from generator? \n");
2302  for(UInt_t imc = 0; imc < 4; imc++)
2303  {
2304  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc )
2305  {
2306  mc->GetCocktailGenerator(clus->GetLabelAt(imc),genName);
2307 
2308  Bool_t generOK = kFALSE;
2309  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
2310  {
2311  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
2312  }
2313 
2314  if ( !generOK )
2315  {
2316  amp-=ampOrg*eDepFrac[imc];
2317 
2318  edepTotFrac-=eDepFrac[imc];
2319 
2320  eDepFrac[imc] = 0;
2321  }
2322 
2323  } // eDep > 0
2324  } // 4 possible loop
2325 
2326  } // accept at least one kind of generator
2327 
2328  //
2329  // Add MC label and Edep to corresponding array (to be used later in digits)
2330  //
2331  Int_t nLabels = 0;
2332  for(UInt_t imc = 0; imc < 4; imc++)
2333  {
2334  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc && edepTotFrac > 0 )
2335  {
2336  nLabels++;
2337 
2338  labeArr.Set(nLabels);
2339  labeArr.AddAt(clus->GetLabelAt(imc), nLabels-1);
2340 
2341  eDepArr.Set(nLabels);
2342  eDepArr.AddAt( (eDepFrac[imc]/edepTotFrac) * amp, nLabels-1);
2343  // use as deposited energy a fraction of the simulated energy (smeared and with noise)
2344  } // edep > 0
2345 
2346  } // mc cell label loop
2347 
2348  //
2349  // If no label found, reject cell
2350  // It can happen to have this case (4 MC labels per cell is not enough for some cases)
2351  // better to remove. To be treated carefully.
2352  //
2353  if ( nLabels == 0 ) amp = 0;
2354 
2355 }
2356 
2364 //______________________________________________________________________________
2365 void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
2366  AliVCaloCells* cells,
2367  AliVCluster* clu)
2368 {
2369  if (!clu)
2370  {
2371  AliInfo("Cluster pointer null!");
2372  return;
2373  }
2374 
2376  else if (fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
2377  else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
2378 }
2379 
2388 //_____________________________________________________________________________________________
2390  AliVCaloCells* cells,
2391  AliVCluster* clu)
2392 {
2393  Double_t eCell = 0.;
2394  Float_t fraction = 1.;
2395  Float_t recalFactor = 1.;
2396 
2397  Int_t absId = -1;
2398  Int_t iTower = -1, iIphi = -1, iIeta = -1;
2399  Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
2400  Float_t weight = 0., totalWeight=0.;
2401  Float_t newPos[3] = {-1.,-1.,-1.};
2402  Double_t pLocal[3], pGlobal[3];
2403  Bool_t shared = kFALSE;
2404 
2405  Float_t clEnergy = clu->E(); //Energy already recalibrated previously
2406  if (clEnergy <= 0) return;
2407 
2408  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
2409  Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
2410 
2411  //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
2412 
2413  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
2414  {
2415  absId = clu->GetCellAbsId(iDig);
2416  fraction = clu->GetCellAmplitudeFraction(iDig);
2417  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2418 
2419  if (!fCellsRecalibrated)
2420  {
2421  geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
2422  geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
2423  if (IsRecalibrationOn()) {
2424  if(fUse1Drecalib)
2425  recalFactor = GetEMCALChannelRecalibrationFactor1D(absId);
2426  else
2427  recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
2428  }
2429  }
2430 
2431  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2432 
2433  weight = GetCellWeight(eCell,clEnergy);
2434  totalWeight += weight;
2435 
2436  geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
2437  //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
2438  geom->GetGlobal(pLocal,pGlobal,iSupModMax);
2439  //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
2440 
2441  for (int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
2442  }// cell loop
2443 
2444  if (totalWeight>0)
2445  {
2446  for (int i=0; i<3; i++ ) newPos[i] /= totalWeight;
2447  }
2448 
2449  //Float_t pos[]={0,0,0};
2450  //clu->GetPosition(pos);
2451  //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
2452  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
2453 
2454  if (iSupModMax > 1) { //sector 1
2455  newPos[0] +=fMisalTransShift[3];//-=3.093;
2456  newPos[1] +=fMisalTransShift[4];//+=6.82;
2457  newPos[2] +=fMisalTransShift[5];//+=1.635;
2458  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
2459  } else { //sector 0
2460  newPos[0] +=fMisalTransShift[0];//+=1.134;
2461  newPos[1] +=fMisalTransShift[1];//+=8.2;
2462  newPos[2] +=fMisalTransShift[2];//+=1.197;
2463  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
2464  }
2465  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
2466 
2467  clu->SetPosition(newPos);
2468 }
2469 
2478 //____________________________________________________________________________________________
2480  AliVCaloCells* cells,
2481  AliVCluster* clu)
2482 {
2483  Double_t eCell = 1.;
2484  Float_t fraction = 1.;
2485  Float_t recalFactor = 1.;
2486 
2487  Int_t absId = -1;
2488  Int_t iTower = -1;
2489  Int_t iIphi = -1, iIeta = -1;
2490  Int_t iSupMod = -1, iSupModMax = -1;
2491  Int_t iphi = -1, ieta =-1;
2492  Bool_t shared = kFALSE;
2493 
2494  Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
2495 
2496  if (clEnergy <= 0)
2497  return;
2498  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
2499  Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
2500 
2501  Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
2502  Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
2503  Int_t startingSM = -1;
2504 
2505  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
2506  {
2507  absId = clu->GetCellAbsId(iDig);
2508  fraction = clu->GetCellAmplitudeFraction(iDig);
2509  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2510 
2511  if (iDig==0) startingSM = iSupMod;
2512  else if (iSupMod != startingSM) areInSameSM = kFALSE;
2513 
2514  eCell = cells->GetCellAmplitude(absId);
2515 
2516  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2517  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
2518 
2519  if (!fCellsRecalibrated)
2520  {
2521  if (IsRecalibrationOn()) {
2522  if(fUse1Drecalib)
2523  recalFactor = GetEMCALChannelRecalibrationFactor1D(absId);
2524  else
2525  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2526  }
2527  }
2528 
2529  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2530 
2531  weight = GetCellWeight(eCell,clEnergy);
2532  if (weight < 0) weight = 0;
2533  totalWeight += weight;
2534  weightedCol += ieta*weight;
2535  weightedRow += iphi*weight;
2536 
2537  //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
2538  }// cell loop
2539 
2540  Float_t xyzNew[]={-1.,-1.,-1.};
2541  if (areInSameSM == kTRUE)
2542  {
2543  //printf("In Same SM\n");
2544  weightedCol = weightedCol/totalWeight;
2545  weightedRow = weightedRow/totalWeight;
2546  geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2547  }
2548  else
2549  {
2550  //printf("In Different SM\n");
2551  geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2552  }
2553 
2554  clu->SetPosition(xyzNew);
2555 }
2556 
2564 //___________________________________________________________________________________________
2566  AliVCaloCells* cells,
2567  AliVCluster * cluster)
2568 {
2569  if (!fRecalDistToBadChannels) return;
2570 
2571  if (!cluster)
2572  {
2573  AliInfo("Cluster pointer null!");
2574  return;
2575  }
2576 
2577  // Get channels map of the supermodule where the cluster is.
2578  Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
2579  Bool_t shared = kFALSE;
2580  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
2581  TH2D* hMap = 0x0;
2582  TH1C* hMap1D = 0x0;
2583 
2584  if(!fUse1Dmap)
2585  hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
2586  else
2587  hMap1D = (TH1C*)fEMCALBadChannelMap->At(0);
2588 
2589  Int_t dRrow, dRcol;
2590  Float_t minDist = 10000.;
2591  Float_t dist = 0.;
2592 
2593  // Loop on tower status map
2594  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2595  {
2596  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2597  {
2598  // Check if tower is bad.
2599  Int_t status=0;
2600  if (fUse1Dmap)
2601  status = hMap1D->GetBinContent(geom->GetAbsCellIdFromCellIndexes(iSupMod, irow, icol));
2602  else
2603  status = hMap->GetBinContent(icol,irow);
2604 
2605  if(status==0) continue;
2606 
2607  //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
2608  // iSupMod,icol, irow, icolM,irowM);
2609 
2610  dRrow=TMath::Abs(irowM-irow);
2611  dRcol=TMath::Abs(icolM-icol);
2612  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2613  if (dist < minDist)
2614  {
2615  //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
2616  minDist = dist;
2617  }
2618  }
2619  }
2620 
2621  // In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
2622  if (shared)
2623  {
2624  TH2D* hMap2 = 0;
2625  TH1C* hMap1D2 = 0;
2626  Int_t iSupMod2 = -1;
2627 
2628  // The only possible combinations are (0,1), (2,3) ... (8,9)
2629  if (iSupMod%2) iSupMod2 = iSupMod-1;
2630  else iSupMod2 = iSupMod+1;
2631 
2632  if(!fUse1Dmap)
2633  hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
2634  else
2635  hMap1D2 = (TH1C*)fEMCALBadChannelMap->At(0);
2636 
2637 
2638  // Loop on tower status map of second super module
2639  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2640  {
2641  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2642  {
2643  // Check if tower is bad.
2644  Int_t status=0;
2645  if (fUse1Dmap)
2646  status = hMap1D2->GetBinContent(geom->GetAbsCellIdFromCellIndexes(iSupMod2, irow, icol));
2647  else
2648  status = hMap2->GetBinContent(icol,irow);
2649 
2650  if(status==0) continue;
2651 
2652  //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",
2653  // iSupMod2,icol, irow,iSupMod,icolM,irowM);
2654  dRrow=TMath::Abs(irow-irowM);
2655 
2656  if (iSupMod%2)
2657  dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
2658  else
2659  dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
2660 
2661  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2662  if (dist < minDist) minDist = dist;
2663  }
2664  }
2665  }// shared cluster in 2 SuperModules
2666 
2667  AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
2668  cluster->SetDistanceToBadChannel(minDist);
2669 }
2670 
2677 //__________________________________________________________________
2678 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
2679 {
2680  if (!cluster)
2681  {
2682  AliInfo("Cluster pointer null!");
2683  return;
2684  }
2685 
2686  if (cluster->GetM02() != 0)
2687  fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
2688 
2689  Float_t pidlist[AliPID::kSPECIESCN+1];
2690  for (Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
2691 
2692  cluster->SetPID(pidlist);
2693 }
2694 
2715 //___________________________________________________________________________________________________________________
2717  AliVCaloCells* cells, AliVCluster * cluster,
2718  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2719  Float_t & enAfterCuts, Float_t & l0, Float_t & l1,
2720  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2721  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2722 {
2723  if (!cluster)
2724  {
2725  AliInfo("Cluster pointer null!");
2726  return;
2727  }
2728 
2729  Double_t eCell = 0.;
2730  Double_t tCell = 0.;
2731  Float_t fraction = 1.;
2732  Float_t recalFactor = 1.;
2733  Bool_t isLowGain = kFALSE;
2734 
2735  Int_t iSupMod = -1;
2736  Int_t iTower = -1;
2737  Int_t iIphi = -1;
2738  Int_t iIeta = -1;
2739  Int_t iphi = -1;
2740  Int_t ieta = -1;
2741  Double_t etai = -1.;
2742  Double_t phii = -1.;
2743 
2744  Int_t nstat = 0 ;
2745  Float_t wtot = 0.;
2746  Double_t w = 0.;
2747  Double_t etaMean = 0.;
2748  Double_t phiMean = 0.;
2749 
2750  Double_t pGlobal[3];
2751 
2752  // Loop on cells, calculate the cluster energy, in case a cut on cell energy is added,
2753  // or the non linearity correction was applied
2754  // and to check if the cluster is between 2 SM in eta
2755  Int_t iSM0 = -1;
2756  Bool_t shared = kFALSE;
2757  Float_t energy = 0;
2758 
2759  enAfterCuts = 0;
2760  l0 = 0; l1 = 0;
2761  disp = 0; dEta = 0; dPhi = 0;
2762  sEta = 0; sPhi = 0; sEtaPhi = 0;
2763 
2764  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2765  {
2766  // Get from the absid the supermodule, tower and eta/phi numbers
2767 
2768  Int_t absId = cluster->GetCellAbsId(iDigit);
2769 
2770  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2771  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2772 
2773  // Check if there are cells of different SM
2774  if (iDigit == 0 ) iSM0 = iSupMod;
2775  else if (iSupMod!= iSM0) shared = kTRUE;
2776 
2777  // Get the cell energy, if recalibration is on, apply factors
2778  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2779  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2780 
2781  if (IsRecalibrationOn()){
2782  if(fUse1Drecalib)
2783  recalFactor = GetEMCALChannelRecalibrationFactor1D(absId);
2784  else
2785  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2786  }
2787 
2788  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2789  tCell = cells->GetCellTime (absId);
2790  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2791 
2792  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2793  tCell*=1e9;
2794  tCell-=fConstantTimeShift;
2795 
2796  if(eCell > 0.05) // at least a minimum 50 MeV cut
2797  energy += eCell;
2798 
2799  if(eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut)
2800  enAfterCuts += eCell;
2801  } // cell loop
2802 
2803 
2804  // Loop on cells to calculate weights and shower shape terms parameters
2805  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2806  {
2807  // Get from the absid the supermodule, tower and eta/phi numbers
2808  Int_t absId = cluster->GetCellAbsId(iDigit);
2809 
2810  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2811  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2812 
2813  //Get the cell energy, if recalibration is on, apply factors
2814  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2815  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2816 
2818  if(fUse1Drecalib)
2819  recalFactor = GetEMCALChannelRecalibrationFactor1D(absId);
2820  else
2821  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2822  }
2823 
2824  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2825  tCell = cells->GetCellTime (absId);
2826  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2827 
2828  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2829  tCell*=1e9;
2830  tCell-=fConstantTimeShift;
2831 
2832  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2833  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2834  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2835 
2836  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2837  {
2838  w = GetCellWeight(eCell, energy);
2839 
2840  // Cell index
2841  if ( fShowerShapeCellLocationType == 0 )
2842  {
2843  etai=(Double_t)ieta;
2844  phii=(Double_t)iphi;
2845  }
2846  // Cell angle location
2847  else if( fShowerShapeCellLocationType == 1 )
2848  {
2849  geom->EtaPhiFromIndex(absId, etai, phii);
2850  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2851  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2852  }
2853  else
2854  {
2855  geom->GetGlobal(absId,pGlobal);
2856 
2857  // Cell x-z location
2858  if( fShowerShapeCellLocationType == 2 )
2859  {
2860  etai = pGlobal[2];
2861  phii = pGlobal[0];
2862  }
2863  // Cell r-z location
2864  else
2865  {
2866  etai = pGlobal[2];
2867  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2868  }
2869  }
2870 
2871  if (w > 0.0)
2872  {
2873  wtot += w ;
2874  nstat++;
2875 
2876  // Shower shape
2877  sEta += w * etai * etai ;
2878  etaMean += w * etai ;
2879  sPhi += w * phii * phii ;
2880  phiMean += w * phii ;
2881  sEtaPhi += w * etai * phii ;
2882  }
2883  }
2884  else if(eCell > 0.05)
2885  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2886  } // cell loop
2887 
2888  //printf("sEta %f sPhi %f etaMean %f phiMean %f sEtaPhi %f wtot %f\n",sEta,sPhi,etaMean,phiMean,sEtaPhi, wtot);
2889 
2890  // Normalize to the weight
2891  if (wtot > 0)
2892  {
2893  etaMean /= wtot ;
2894  phiMean /= wtot ;
2895  }
2896  else
2897  AliDebug(2,Form("Wrong weight %f\n", wtot));
2898 
2899  // Loop on cells to calculate dispersion
2900  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2901  {
2902  // Get from the absid the supermodule, tower and eta/phi numbers
2903  Int_t absId = cluster->GetCellAbsId(iDigit);
2904 
2905  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2906  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2907 
2908  //Get the cell energy, if recalibration is on, apply factors
2909  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2910  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2911  if (IsRecalibrationOn()) {
2912  if(fUse1Drecalib)
2913  recalFactor = GetEMCALChannelRecalibrationFactor1D(absId);
2914  else
2915  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2916  }
2917 
2918  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2919  tCell = cells->GetCellTime (absId);
2920  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2921 
2922  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2923  tCell*=1e9;
2924  tCell-=fConstantTimeShift;
2925 
2926  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2927  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2928  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2929 
2930  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2931  {
2932  w = GetCellWeight(eCell,cluster->E());
2933 
2934  // Cell index
2935  if ( fShowerShapeCellLocationType == 0 )
2936  {
2937  etai=(Double_t)ieta;
2938  phii=(Double_t)iphi;
2939  }
2940  // Cell angle location
2941  else if( fShowerShapeCellLocationType == 1 )
2942  {
2943  geom->EtaPhiFromIndex(absId, etai, phii);
2944  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2945  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2946  }
2947  else
2948  {
2949  geom->GetGlobal(absId,pGlobal);
2950 
2951  // Cell x-z location
2952  if( fShowerShapeCellLocationType == 2 )
2953  {
2954  etai = pGlobal[2];
2955  phii = pGlobal[0];
2956  }
2957  // Cell r-z location
2958  else
2959  {
2960  etai = pGlobal[2];
2961  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2962  }
2963  }
2964 
2965  if (w > 0.0)
2966  {
2967  disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
2968  dEta += w * (etai-etaMean)*(etai-etaMean) ;
2969  dPhi += w * (phii-phiMean)*(phii-phiMean) ;
2970  }
2971  }
2972  else
2973  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2974  }// cell loop
2975 
2976  // Normalize to the weigth and set shower shape parameters
2977  if (wtot > 0 && nstat > 1)
2978  {
2979  disp /= wtot ;
2980  dEta /= wtot ;
2981  dPhi /= wtot ;
2982  sEta /= wtot ;
2983  sPhi /= wtot ;
2984  sEtaPhi /= wtot ;
2985 
2986  sEta -= etaMean * etaMean ;
2987  sPhi -= phiMean * phiMean ;
2988  sEtaPhi -= etaMean * phiMean ;
2989 
2990  l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2991  l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2992 
2993  //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);
2994 
2995  }
2996  else
2997  {
2998  l0 = 0. ;
2999  l1 = 0. ;
3000  dEta = 0. ; dPhi = 0. ; disp = 0. ;
3001  sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
3002  }
3003 }
3004 
3023 //___________________________________________________________________________________________________________________
3025  AliVCaloCells* cells, AliVCluster * cluster,
3026  Float_t & l0, Float_t & l1,
3027  Float_t & disp, Float_t & dEta, Float_t & dPhi,
3028  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
3029 {
3030  Float_t newEnergy = 0;
3031  Float_t cellEmin = 0.05; // 50 MeV
3032  Float_t cellTimeWindow = 1000; // open cut
3033  Int_t bc = 0;
3035  cellEmin, cellTimeWindow, bc,
3036  newEnergy, l0, l1, disp,
3037  dEta, dPhi, sEta, sPhi, sEtaPhi);
3038 }
3039 
3048 //____________________________________________________________________________________________
3050  AliVCaloCells* cells,
3051  AliVCluster * cluster)
3052 {
3053  Float_t l0 = 0., l1 = 0.;
3054  Float_t disp = 0., dEta = 0., dPhi = 0.;
3055  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
3056 
3057  AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
3058  dEta, dPhi, sEta, sPhi, sEtaPhi);
3059 
3060  cluster->SetM02(l0);
3061  cluster->SetM20(l1);
3062  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
3063 
3064 }
3065 
3078 //____________________________________________________________________________________________
3080  AliVCaloCells* cells, AliVCluster * cluster,
3081  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
3082  Float_t & enAfterCuts)
3083 {
3084  Float_t l0 = 0., l1 = 0.;
3085  Float_t disp = 0., dEta = 0., dPhi = 0.;
3086  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
3087 
3089  cellEcut, cellTimeCut, bc,
3090  enAfterCuts, l0, l1, disp,
3091  dEta, dPhi, sEta, sPhi, sEtaPhi);
3092 
3093  cluster->SetM02(l0);
3094  cluster->SetM20(l1);
3095  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
3096 
3097 }
3098 
3112 //____________________________________________________________________________
3113 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
3114  TObjArray * clusterArr,
3115  const AliEMCALGeometry *geom,
3116  AliMCEvent * mc)
3117 {
3118  fMatchedTrackIndex ->Reset();
3119  fMatchedClusterIndex->Reset();
3120  fResidualPhi->Reset();
3121  fResidualEta->Reset();
3122 
3123  fMatchedTrackIndex ->Set(1000);
3124  fMatchedClusterIndex->Set(1000);
3125  fResidualPhi->Set(1000);
3126  fResidualEta->Set(1000);
3127 
3128  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
3129  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
3130 
3131  // Init the magnetic field if not already on
3132  if (!TGeoGlobalMagField::Instance()->GetField())
3133  {
3134  if (!event->InitMagneticField())
3135  {
3136  AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
3137  }
3138  } // Init mag field
3139 
3140  if (esdevent)
3141  {
3142  UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
3143  UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
3144  Bool_t desc1 = (mask1 >> 3) & 0x1;
3145  Bool_t desc2 = (mask2 >> 3) & 0x1;
3146  if (desc1==0 || desc2==0) {
3147  // AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
3148  // mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
3149  // mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
3150  fITSTrackSA=kTRUE;
3151  }
3152  }
3153 
3154  TObjArray *clusterArray = 0x0;
3155  if (!clusterArr)
3156  {
3157  clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
3158  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
3159  {
3160  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
3161  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells()))
3162  continue;
3163  clusterArray->AddAt(cluster,icl);
3164  }
3165  }
3166 
3167  Int_t matched=0;
3168  Double_t cv[21];
3169  TString genName;
3170  for (Int_t i=0; i<21;i++) cv[i]=0;
3171  for (Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
3172  {
3173  AliExternalTrackParam *trackParam = 0;
3174  Int_t mcLabel = -1;
3175  // If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
3176  AliESDtrack *esdTrack = 0;
3177  AliAODTrack *aodTrack = 0;
3178  if (esdevent)
3179  {
3180  esdTrack = esdevent->GetTrack(itr);
3181  if (!esdTrack) continue;
3182  if (!IsAccepted(esdTrack)) continue;
3183  if (esdTrack->Pt()<fCutMinTrackPt) continue;
3184 
3185  if ( TMath::Abs(esdTrack->Eta()) > 0.9 ) continue;
3186 
3187  // Save some time and memory in case of no DCal present
3188  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
3189  {
3190  Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
3191  if ( phi <= 10 || phi >= 250 ) continue;
3192  }
3193 
3194  if (!fITSTrackSA) // if TPC Available
3195  {
3196  if ( fUseOuterTrackParam )
3197  trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetOuterParam());
3198  else
3199  trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
3200  }
3201  else
3202  trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
3203 
3204  mcLabel = TMath::Abs(esdTrack->GetLabel());
3205  }
3206 
3207  // If the input event is AOD, the starting point for extrapolation is at vertex
3208  // AOD tracks are selected according to its filterbit.
3209  else if (aodevent)
3210  {
3211  aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(itr));
3212 
3213  if (!aodTrack) AliFatal("Not a standard AOD");
3214 
3215  if (!aodTrack) continue;
3216 
3217  if (fAODTPCOnlyTracks)
3218  { // Match with TPC only tracks, default from May 2013, before filter bit 32
3219  //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
3220  if (!aodTrack->IsTPCConstrained()) continue ;
3221  }
3222  else if (fAODHybridTracks)
3223  { // Match with hybrid tracks
3224  //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
3225  if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
3226  } else
3227  { // Match with tracks on a mask
3228  //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask));
3229  if (!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks
3230  }
3231 
3232  if (aodTrack->Pt() < fCutMinTrackPt) continue;
3233 
3234  if ( TMath::Abs(aodTrack->Eta()) > 0.9 ) continue;
3235 
3236  // Save some time and memory in case of no DCal present
3237  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
3238  {
3239  Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
3240  if ( phi <= 10 || phi >= 250 ) continue;
3241  }
3242 
3243  Double_t pos[3],mom[3];
3244 
3245  if ( fUseTrackDCA )
3246  aodTrack->GetXYZ(pos);
3247  else
3248  aodTrack->XvYvZv(pos);
3249 
3250  aodTrack->GetPxPyPz(mom);
3251  AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",
3252  itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
3253 
3254  trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
3255 
3256  mcLabel = TMath::Abs(aodTrack->GetLabel());
3257  }
3258 
3259  //Return if the input data is not "AOD" or "ESD"
3260  else
3261  {
3262  AliWarning("Wrong input data type! Should be \"AOD\" or \"ESD\" ");
3263  if (clusterArray)
3264  {
3265  clusterArray->Clear();
3266  delete clusterArray;
3267  }
3268  return;
3269  }
3270 
3271  if (!trackParam) continue;
3272 
3273  //
3274  // Check if track comes from a particular MC generator, do not include it
3275  // if it is not a selected one
3276  //
3277  if( mc && fMCGenerToAcceptForTrack && fNMCGenerToAccept > 0 )
3278  {
3279  mc->GetCocktailGenerator(mcLabel,genName);
3280 
3281  Bool_t generOK = kFALSE;
3282  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
3283  {
3284  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
3285  }
3286 
3287  if ( !generOK ) continue;
3288  }
3289 
3290  // Extrapolate the track to EMCal surface, see AliEMCALRecoUtilsBase
3291  AliExternalTrackParam emcalParam(*trackParam);
3292  Float_t eta, phi, pt;
3293  if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
3294  {
3295  if (aodevent && trackParam) delete trackParam;
3296  if (fITSTrackSA && trackParam) delete trackParam;
3297  continue;
3298  }
3299 
3300  if ( TMath::Abs(eta) > 0.75 )
3301  {
3302  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
3303  continue;
3304  }
3305 
3306  // Save some time and memory in case of no DCal present
3307  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
3308  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
3309  {
3310  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
3311  continue;
3312  }
3313 
3314  //Find matched clusters
3315  Int_t index = -1;
3316  Float_t dEta = -999, dPhi = -999;
3317  if (!clusterArr)
3318  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
3319  else
3320  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr , dEta, dPhi);
3321 
3322 
3323  if (index>-1)
3324  {
3325  fMatchedTrackIndex ->AddAt(itr,matched);
3326  fMatchedClusterIndex ->AddAt(index,matched);
3327  fResidualEta ->AddAt(dEta,matched);
3328  fResidualPhi ->AddAt(dPhi,matched);
3329  matched++;
3330  }
3331 
3332  if (aodevent && trackParam) delete trackParam;
3333  if (fITSTrackSA && trackParam) delete trackParam;
3334  }//track loop
3335 
3336  if (clusterArray)
3337  {
3338  clusterArray->Clear();
3339  delete clusterArray;
3340  }
3341 
3342  AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
3343 
3344  fMatchedTrackIndex ->Set(matched);
3345  fMatchedClusterIndex ->Set(matched);
3346  fResidualPhi ->Set(matched);
3347  fResidualEta ->Set(matched);
3348 }
3349 
3361 //________________________________________________________________________________
3363  const AliVEvent *event,
3364  const AliEMCALGeometry *geom,
3365  Float_t &dEta, Float_t &dPhi)
3366 {
3367  Int_t index = -1;
3368 
3369  if ( TMath::Abs(track->Eta()) > 0.9 ) return index;
3370 
3371  // Save some time and memory in case of no DCal present
3372  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
3373  {
3374  Double_t phiV = track->Phi()*TMath::RadToDeg();
3375  if ( phiV <= 10 || phiV >= 250 ) return index;
3376  }
3377 
3378  AliExternalTrackParam *trackParam = 0;
3379  if (!fITSTrackSA) // If TPC
3380  {
3381  if ( fUseOuterTrackParam )
3382  trackParam = const_cast<AliExternalTrackParam*>(track->GetOuterParam());
3383  else
3384  trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
3385  }
3386  else
3387  trackParam = new AliExternalTrackParam(*track);
3388 
3389  if (!trackParam) return index;
3390  AliExternalTrackParam emcalParam(*trackParam);
3391 
3392  Float_t eta, phi, pt;
3393  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
3394  {
3395  if (fITSTrackSA) delete trackParam;
3396  return index;
3397  }
3398 
3399  if ( TMath::Abs(eta) > 0.75 )
3400  {
3401  if (fITSTrackSA) delete trackParam;
3402  return index;
3403  }
3404 
3405  // Save some time and memory in case of no DCal present
3406  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
3407  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
3408  {
3409  if (fITSTrackSA) delete trackParam;
3410  return index;
3411  }
3412 
3413  TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
3414 
3415  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
3416  {
3417  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
3418  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
3419  clusterArr->AddAt(cluster,icl);
3420  }
3421 
3422  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
3423  clusterArr->Clear();
3424  delete clusterArr;
3425  if (fITSTrackSA) delete trackParam;
3426 
3427  return index;
3428 }
3429 
3440 //_______________________________________________________________________________________________
3441 Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
3442  AliExternalTrackParam *trkParam,
3443  const TObjArray * clusterArr,
3444  Float_t &dEta, Float_t &dPhi)
3445 {
3446  dEta=-999, dPhi=-999;
3447  Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
3448  Int_t index = -1;
3449  Float_t tmpEta=-999, tmpPhi=-999;
3450 
3451  Double_t exPos[3] = {0.,0.,0.};
3452  if (!emcalParam->GetXYZ(exPos)) return index;
3453 
3454  Float_t clsPos[3] = {0.,0.,0.};
3455  for (Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
3456  {
3457  AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
3458 
3459  if (!cluster || !cluster->IsEMCAL()) continue;
3460 
3461  cluster->GetPosition(clsPos);
3462 
3463  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));
3464  if (dR > fClusterWindow) continue;
3465 
3466  AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
3467 
3468  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
3469 
3470  if (fCutEtaPhiSum)
3471  {
3472  Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
3473  if (tmpR<dRMax)
3474  {
3475  dRMax=tmpR;
3476  dEtaMax=tmpEta;
3477  dPhiMax=tmpPhi;
3478  index=icl;
3479  }
3480  }
3481  else if (fCutEtaPhiSeparate)
3482  {
3483  if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
3484  {
3485  dEtaMax = tmpEta;
3486  dPhiMax = tmpPhi;
3487  index=icl;
3488  }
3489  }
3490  else
3491  {
3492  AliError("Please specify your cut criteria");
3493  AliError("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()");
3494  AliError("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()");
3495  return index;
3496  }
3497  }
3498 
3499  dEta=dEtaMax;
3500  dPhi=dPhiMax;
3501 
3502  return index;
3503 }
3504 
3517 //---------------------------------------------------------------------------------
3518 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
3519  const AliVCluster *cluster,
3520  Float_t &tmpEta,
3521  Float_t &tmpPhi)
3522 {
3523  return AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
3524 }
3525 
3535 //_______________________________________________________________________
3537  Float_t &dEta, Float_t &dPhi)
3538 {
3539  if (FindMatchedPosForCluster(clsIndex) >= 999)
3540  {
3541  AliDebug(2,"No matched tracks found!\n");
3542  dEta=999.;
3543  dPhi=999.;
3544  return;
3545  }
3546 
3547  dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
3548  dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
3549 }
3550 
3560 //______________________________________________________________________________________________
3562 {
3563  if (FindMatchedPosForTrack(trkIndex) >= 999)
3564  {
3565  AliDebug(2,"No matched cluster found!\n");
3566  dEta=999.;
3567  dPhi=999.;
3568  return;
3569  }
3570 
3571  dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
3572  dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
3573 }
3574 
3581 //__________________________________________________________
3583 {
3584  if (IsClusterMatched(clsIndex))
3585  return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
3586  else
3587  return -1;
3588 }
3589 
3596 //__________________________________________________________
3598 {
3599  if (IsTrackMatched(trkIndex))
3600  return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
3601  else
3602  return -1;
3603 }
3604 
3612 //______________________________________________________________
3614 {
3615  if (FindMatchedPosForCluster(clsIndex) < 999)
3616  return kTRUE;
3617  else
3618  return kFALSE;
3619 }
3620 
3628 //____________________________________________________________
3630 {
3631  if (FindMatchedPosForTrack(trkIndex) < 999)
3632  return kTRUE;
3633  else
3634  return kFALSE;
3635 }
3636 
3644 //______________________________________________________________________
3646 {
3647  Float_t tmpR = fCutR;
3648  UInt_t pos = 999;
3649 
3650  for (Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
3651  {
3652  if (fMatchedClusterIndex->At(i)==clsIndex)
3653  {
3654  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3655  if (r<tmpR)
3656  {
3657  pos=i;
3658  tmpR=r;
3659  AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3660  fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3661  }
3662  }
3663  }
3664  return pos;
3665 }
3666 
3674 //____________________________________________________________________
3676 {
3677  Float_t tmpR = fCutR;
3678  UInt_t pos = 999;
3679 
3680  for (Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
3681  {
3682  if (fMatchedTrackIndex->At(i)==trkIndex)
3683  {
3684  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3685  if (r<tmpR)
3686  {
3687  pos=i;
3688  tmpR=r;
3689  AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3690  fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3691  }
3692  }
3693  }
3694  return pos;
3695 }
3696 
3711 //__________________________________________________________________________
3713  const AliEMCALGeometry *geom,
3714  AliVCaloCells* cells, Int_t bc)
3715 {
3716  Bool_t isGood=kTRUE;
3717 
3718  if (!cluster || !cluster->IsEMCAL()) return kFALSE;
3719  if (ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
3720  if (!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
3721  if (IsExoticCluster(cluster, cells,bc)) return kFALSE;
3722 
3723  return isGood;
3724 }
3725 
3736 //__________________________________________________________
3738 {
3739  UInt_t status = esdTrack->GetStatus();
3740 
3741  Int_t nClustersITS = esdTrack->GetITSclusters(0);
3742  Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
3743 
3744  Float_t chi2PerClusterITS = -1;
3745  Float_t chi2PerClusterTPC = -1;
3746  if (nClustersITS!=0)
3747  chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
3748  if (nClustersTPC!=0)
3749  chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
3750 
3751  //
3752  // DCA cuts
3753  // Only to be used for primary
3754  //
3755  if ( fTrackCutsType == kGlobalCut )
3756  {
3757  Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01);
3758  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
3759 
3760  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3761 
3762  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3763  }
3764  else if( fTrackCutsType == kGlobalCut2011 )
3765  {
3766  Float_t maxDCAToVertexXYPtDep = 0.0105 + 0.0350/TMath::Power(esdTrack->Pt(),1.1);
3767  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2011()
3768 
3769  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3770 
3771  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3772  }
3773 
3774  Float_t b[2];
3775  Float_t bCov[3];
3776  esdTrack->GetImpactParameters(b,bCov);
3777  if (bCov[0]<=0 || bCov[2]<=0)
3778  {
3779  AliDebug(1, "Estimated b resolution lower or equal zero!");
3780  bCov[0]=0; bCov[2]=0;
3781  }
3782 
3783  Float_t dcaToVertexXY = b[0];
3784  Float_t dcaToVertexZ = b[1];
3785  Float_t dcaToVertex = -1;
3786 
3787  if (fCutDCAToVertex2D)
3788  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY / fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY +
3789  dcaToVertexZ*dcaToVertexZ / fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ );
3790  else
3791  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
3792 
3793  // cut the track?
3794 
3795  Bool_t cuts[kNCuts];
3796  for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
3797 
3798  // track quality cuts
3799  if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
3800  cuts[0]=kTRUE;
3801  if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
3802  cuts[1]=kTRUE;
3803  if (nClustersTPC<fCutMinNClusterTPC)
3804  cuts[2]=kTRUE;
3805  if (nClustersITS<fCutMinNClusterITS)
3806  cuts[3]=kTRUE;
3807  if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
3808  cuts[4]=kTRUE;
3809  if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
3810  cuts[5]=kTRUE;
3811  if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
3812  cuts[6]=kTRUE;
3813  if (fCutDCAToVertex2D && dcaToVertex > 1)
3814  cuts[7] = kTRUE;
3815  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
3816  cuts[8] = kTRUE;
3817  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
3818  cuts[9] = kTRUE;
3819 
3821  {
3822  //Require at least one SPD point + anything else in ITS
3823  if ( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
3824  cuts[10] = kTRUE;
3825  }
3826 
3827  // ITS
3829  {
3830  if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin))
3831  {
3832  // TPC tracks
3833  cuts[11] = kTRUE;
3834  }
3835  else
3836  {
3837  // ITS standalone tracks
3839  {
3840  if (status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
3841  }
3842  else if (fCutRequireITSpureSA)
3843  {
3844  if (!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
3845  }
3846  }
3847  }
3848 
3849  Bool_t cut=kFALSE;
3850  for (Int_t i=0; i<kNCuts; i++)
3851  if (cuts[i]) { cut = kTRUE ; }
3852 
3853  // cut the track
3854  if (cut)
3855  return kFALSE;
3856  else
3857  return kTRUE;
3858 }
3859 
3865 //_____________________________________
3867 {
3868  switch (fTrackCutsType)
3869  {
3870  case kTPCOnlyCut:
3871  {
3872  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()"));
3873  //TPC
3874  SetMinNClustersTPC(70);
3876  SetAcceptKinkDaughters(kFALSE);
3877  SetRequireTPCRefit(kFALSE);
3878 
3879  //ITS
3880  SetRequireITSRefit(kFALSE);
3881  SetMaxDCAToVertexZ(3.2);
3882  SetMaxDCAToVertexXY(2.4);
3883  SetDCAToVertex2D(kTRUE);
3884 
3885  break;
3886  }
3887 
3888  case kGlobalCut:
3889  {
3890  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE)"));
3891  //TPC
3892  SetMinNClustersTPC(70);
3894  SetAcceptKinkDaughters(kFALSE);
3895  SetRequireTPCRefit(kTRUE);
3896 
3897  //ITS
3898  SetRequireITSRefit(kTRUE);
3899  SetMaxDCAToVertexZ(2);
3901  SetDCAToVertex2D(kFALSE);
3902 
3903  break;
3904  }
3905 
3906  case kLooseCut:
3907  {
3908  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
3909  SetMinNClustersTPC(50);
3910  SetAcceptKinkDaughters(kTRUE);
3911 
3912  break;
3913  }
3914 
3915  case kITSStandAlone:
3916  {
3917  AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
3918  SetRequireITSRefit(kTRUE);
3919  SetRequireITSStandAlone(kTRUE);
3920  SetITSTrackSA(kTRUE);
3921  break;
3922  }
3923 
3924  case kGlobalCut2011:
3925  {
3926  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE)"));
3927  //TPC
3928  SetMinNClustersTPC(50);
3930  SetAcceptKinkDaughters(kFALSE);
3931  SetRequireTPCRefit(kTRUE);
3932 
3933  //ITS
3934  SetRequireITSRefit(kTRUE);
3935  SetMaxDCAToVertexZ(2);
3937  SetDCAToVertex2D(kFALSE);
3938 
3939  break;
3940  }
3941 
3942  case kLooseCutWithITSrefit:
3943  {
3944  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut plus ITSrefit"));
3945  SetMinNClustersTPC(50);
3946  SetAcceptKinkDaughters(kTRUE);
3947  SetRequireITSRefit(kTRUE);
3948 
3949  break;
3950  }
3951  }
3952 }
3953 
3959 //________________________________________________________________________
3961 {
3962  Int_t nTracks = event->GetNumberOfTracks();
3963  for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
3964  {
3965  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
3966  if (!track)
3967  {
3968  AliWarning(Form("Could not receive track %d", iTrack));
3969  continue;
3970  }
3971 
3972  Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
3973  track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
3974 
3975  /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
3976  AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
3977  if (esdtrack)
3978  {
3979  if (matchClusIndex != -1)
3980  esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
3981  else
3982  esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3983  }
3984  else
3985  {
3986  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
3987  if (matchClusIndex != -1)
3988  aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
3989  else
3990  aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3991  }
3992  }
3993  AliDebug(2,"Track matched to closest cluster");
3994 }
3995 
4002 //_________________________________________________________________________
4004 {
4005  for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
4006  {
4007  AliVCluster *cluster = event->GetCaloCluster(iClus);
4008  if (!cluster->IsEMCAL())
4009  continue;
4010 
4011  //
4012  // Remove old matches in cluster
4013  //
4014  if(cluster->GetNTracksMatched() > 0)
4015  {
4016  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
4017  {
4018  TArrayI arrayTrackMatched(0);
4019  ((AliESDCaloCluster*)cluster)->AddTracksMatched(arrayTrackMatched);
4020  }
4021  else
4022  {
4023  for(Int_t iTrack = 0; iTrack < cluster->GetNTracksMatched(); iTrack++)
4024  {
4025  ((AliAODCaloCluster*)cluster)->RemoveTrackMatched((TObject*)((AliAODCaloCluster*)cluster)->GetTrackMatched(iTrack));
4026  }
4027  }
4028  }
4029 
4030  //
4031  // Find new matches and put them in the cluster
4032  //
4033  Int_t nTracks = event->GetNumberOfTracks();
4034  TArrayI arrayTrackMatched(nTracks);
4035 
4036  // Get the closest track matched to the cluster
4037  Int_t nMatched = 0;
4038  Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
4039  if (matchTrackIndex != -1)
4040  {
4041  arrayTrackMatched[nMatched] = matchTrackIndex;
4042  nMatched++;
4043  }
4044 
4045  // Get all other tracks matched to the cluster
4046  for (Int_t iTrk=0; iTrk<nTracks; ++iTrk)
4047  {
4048  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
4049 
4050  if( !track ) continue;
4051 
4052  if ( iTrk == matchTrackIndex ) continue;
4053 
4054  if ( track->GetEMCALcluster() == iClus )
4055  {
4056  arrayTrackMatched[nMatched] = iTrk;
4057  ++nMatched;
4058  }
4059  }
4060 
4061  AliDebug(2,Form("cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]));
4062 
4063  arrayTrackMatched.Set(nMatched);
4064  AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
4065  if (esdcluster)
4066  esdcluster->AddTracksMatched(arrayTrackMatched);
4067  else if ( nMatched > 0 )
4068  {
4069  AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
4070  if (aodcluster)
4071  {
4072  aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
4073  //AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(event->GetTrack(arrayTrackMatched.At(0)));
4074  //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());
4075  //printf("With specs: pt %.4f, eta %.4f, phi %.4f\n",aodtrack->Pt(),aodtrack->Eta(), aodtrack->Phi());
4076  }
4077  }
4078 
4079  Float_t eta= -999, phi = -999;
4080  if (matchTrackIndex != -1)
4081  GetMatchedResiduals(iClus, eta, phi);
4082 
4083  cluster->SetTrackDistance(phi, eta);
4084  }
4085 
4086  AliDebug(2,"Cluster matched to tracks");
4087 }
4088 
4091  else {
4092  fEMCALRecalibrationFactors = new TObjArray(map->GetEntries());
4093  fEMCALRecalibrationFactors->SetOwner(true);
4094  }
4095  if(!fEMCALRecalibrationFactors->IsOwner()){
4096  // Must claim ownership since the new objects are owend by this instance
4097  fEMCALRecalibrationFactors->SetOwner(kTRUE);
4098  }
4099 
4100  if(!fUse1Drecalib){
4101  for(int i = 0; i < map->GetEntries(); i++){
4102  TH2F *hist = dynamic_cast<TH2F *>(map->At(i));
4103  if(!hist) continue;
4104  this->SetEMCALChannelRecalibrationFactors(i, hist);
4105  }
4106  }else{
4107  TH1S *hist = dynamic_cast<TH1S *>(map->At(0));
4109  }
4110 
4111 }
4112 
4116  fEMCALRecalibrationFactors->SetOwner(true);
4117  }
4118  if(fEMCALRecalibrationFactors->GetEntries() <= iSM) fEMCALRecalibrationFactors->Expand(iSM+1);
4119  if(fEMCALRecalibrationFactors->At(iSM)) fEMCALRecalibrationFactors->RemoveAt(iSM);
4120  TH2F *clone = new TH2F(*h);
4121  clone->SetDirectory(NULL);
4122  fEMCALRecalibrationFactors->AddAt(clone,iSM);
4123 }
4124 
4128  fEMCALRecalibrationFactors->SetOwner(true);
4129  }
4131  TH1S *clone = new TH1S(*h);
4132  clone->SetDirectory(NULL);
4133  fEMCALRecalibrationFactors->AddAt(clone,0);
4134 }
4135 
4137 
4138  if(!fUse1Drecalib){
4139  return (TH2F*)fEMCALRecalibrationFactors->At(iSM);
4140  }else{
4141  const Double_t lowerLimit[]={0,1152,2304,3456,4608,5760,6912,8064,9216,10368,11520,11904,12288,13056,13824,14592,15360,16128,16896,17280};
4142  const Double_t upperLimit[]={1151 ,2303 ,3455 ,4607 ,5759 ,6911 ,8063 ,9215 ,10367,11519,11903,12287,13055,13823,14591,15359,16127,16895,17279,17663};
4143 
4144  TH2F* hist = new TH2F(Form("EMCALRecalFactors_SM%d",iSM),Form("EMCALRecalFactors_SM%d",iSM), 48, 0, 48, 24, 0, 24);
4145  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
4146  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
4147 
4148  for(Int_t iCell=lowerLimit[iSM]; iCell<upperLimit[iSM]; iCell++){
4149  geom->GetCellIndex(iCell,imod,iTower,iIphi,iIeta);
4150  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
4151  hist->SetBinContent(ieta, iphi, ((TH1S*)fEMCALRecalibrationFactors->At(0))->GetBinContent(iCell)/1000);
4152  }
4153  return hist;
4154  }
4155 
4156 }
4157 
4160  else {
4161  fEMCALBadChannelMap = new TObjArray(map->GetEntries());
4162  fEMCALBadChannelMap->SetOwner(true);
4163  }
4164  if(!fEMCALBadChannelMap->IsOwner()) {
4165  // Must claim ownership since the new objects are owend by this instance
4166  fEMCALBadChannelMap->SetOwner(true);
4167  }
4168 
4169  if(!fUse1Dmap){
4170  for(int i = 0; i < map->GetEntries(); i++){
4171  TH2I *hist = dynamic_cast<TH2I *>(map->At(i));
4172  if(!hist) continue;
4173  this->SetEMCALChannelStatusMap(i, hist);
4174  }
4175  }else{
4176  TH1C *hist = dynamic_cast<TH1C *>(map->At(0));
4177  this->SetEMCALChannelStatusMap1D(hist);
4178  }
4179 
4180 }
4181 
4183  if(!fEMCALBadChannelMap){
4184  fEMCALBadChannelMap = new TObjArray(iSM);
4185  fEMCALBadChannelMap->SetOwner(true);
4186  }
4187  if(fEMCALBadChannelMap->GetEntries() <= iSM) fEMCALBadChannelMap->Expand(iSM+1);
4188  if(fEMCALBadChannelMap->At(iSM)) fEMCALBadChannelMap->RemoveAt(iSM);
4189  TH2I *clone = new TH2I(*h);
4190  clone->SetDirectory(NULL);
4191  fEMCALBadChannelMap->AddAt(clone,iSM);
4192 }
4193 
4195  fUse1Dmap = kTRUE;
4196  if(!fEMCALBadChannelMap){
4197  fEMCALBadChannelMap = new TObjArray(1);
4198  fEMCALBadChannelMap->SetOwner(true);
4199  }
4200  if(fEMCALBadChannelMap->At(0)) fEMCALBadChannelMap->RemoveAt(0);
4201  TH1C *clone = new TH1C(*h);
4202  clone->SetDirectory(NULL);
4203  fEMCALBadChannelMap->AddAt(clone,0);
4204 }
4205 
4207 
4208  if(!fUse1Dmap){
4209  return (TH2I*)fEMCALBadChannelMap->At(iSM) ;
4210  }else{
4211  const Double_t lowerLimit[]={0,1152,2304,3456,4608,5760,6912,8064,9216,10368,11520,11904,12288,13056,13824,14592,15360,16128,16896,17280};
4212  const Double_t upperLimit[]={1151 ,2303 ,3455 ,4607 ,5759 ,6911 ,8063 ,9215 ,10367,11519,11903,12287,13055,13823,14591,15359,16127,16895,17279,17663};
4213 
4214  TH2I* hist = new TH2I(Form("EMCALBadChannelMap_Mod%d",iSM),Form("EMCALBadChannelMap_Mod%d",iSM), 48, 0, 48, 24, 0, 24);
4215  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
4216  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
4217 
4218  for(Int_t iCell=lowerLimit[iSM]; iCell<upperLimit[iSM]; iCell++){
4219  geom->GetCellIndex(iCell,imod,iTower,iIphi,iIeta);
4220  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
4221  hist->SetBinContent(ieta, iphi, (Int_t)((TH1C*)fEMCALBadChannelMap->At(0))->GetBinContent(iCell));
4222  }
4223  return hist;
4224  }
4225 }
4226 
4229  else {
4230  fEMCALTimeRecalibrationFactors = new TObjArray(map->GetEntries());
4231  fEMCALTimeRecalibrationFactors->SetOwner(true);
4232  }
4233  if(!fEMCALTimeRecalibrationFactors->IsOwner()) {
4234  // Must claim ownership since the new objects are owend by this instance
4235  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
4236  }
4237  for(int i = 0; i < map->GetEntries(); i++){
4238  TH1F *hist = dynamic_cast<TH1F *>(map->At(i));
4239  if(!hist) continue;
4241  }
4242 }
4243 
4247  fEMCALTimeRecalibrationFactors->SetOwner(true);
4248  }
4249  if(fEMCALTimeRecalibrationFactors->GetEntries() <= bc) fEMCALTimeRecalibrationFactors->Expand(bc+1);
4251  if(fDoUseMergedBC){
4252  TH1S *clone = new TH1S(*(TH1S*)h);
4253  clone->SetDirectory(NULL);
4254  fEMCALTimeRecalibrationFactors->AddAt(clone,bc);
4255  }else{
4256  TH1F *clone = new TH1F(*(TH1F*)h);
4257  clone->SetDirectory(NULL);
4258  fEMCALTimeRecalibrationFactors->AddAt(clone,bc);
4259  }
4260 }
4261 
4264  else {
4265  fEMCALL1PhaseInTimeRecalibration = new TObjArray(map->GetEntries());
4266  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
4267  }
4268  if(!fEMCALL1PhaseInTimeRecalibration->IsOwner()) {
4269  // Must claim ownership since the new objects are owend by this instance
4270  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
4271  }
4272  for(int i = 0; i < map->GetEntries(); i++) {
4273  TH1C *hist = dynamic_cast<TH1C *>(map->At(i));
4274  if(!hist) continue;
4276  }
4277 }
4278 
4282  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
4283  }
4284  if(fEMCALL1PhaseInTimeRecalibration->GetEntries()<parNumber+1) fEMCALL1PhaseInTimeRecalibration->Expand(parNumber+1);
4285  TH1C *clone = new TH1C(*h);
4286  clone->SetDirectory(NULL);
4287  fEMCALL1PhaseInTimeRecalibration->AddAt(clone,parNumber);
4288 }
4289 
4293 //___________________________________________________
4295 {
4296  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
4297  printf("AliEMCALRecoUtils Settings: \n");
4298  printf("\tMisalignment shifts\n");
4299  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,
4301  fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
4302  printf("\tNon linearity function %d, parameters:\n", fNonLinearityFunction);
4303  if (fNonLinearityFunction != 3) // print only if not kNoCorrection
4304  for (Int_t i=0; i<10; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
4305 
4306  printf("\tPosition Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
4307 
4308  printf("\tMatching criteria: ");
4309  if (fCutEtaPhiSum) {
4310  printf("\t\tsqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
4311  } else if (fCutEtaPhiSeparate) {
4312  printf("\t\tdEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
4313  } else {
4314  printf("\t\tError\n");
4315  printf("\t\tplease specify your cut criteria\n");
4316  printf("\t\tTo cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
4317  printf("\t\tTo cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
4318  }
4319 
4320  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);
4321  printf("\tCluster selection window: dR < %2.0f\n",fClusterWindow);
4322 
4323  printf("\tTrack cuts: \n");
4324  printf("\t\tMinimum track pT: %1.2f\n",fCutMinTrackPt);
4325  printf("\t\tAOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask);
4326  printf("\t\tTPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
4327  printf("\t\tAcceptKinks = %d\n",fCutAcceptKinkDaughters);
4328  printf("\t\tMinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
4329  printf("\t\tMaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
4330  printf("\t\tDCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
4331  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
4332 }
void SetEMCALChannelStatusMap1D(const TH1C *h)
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.
Float_t GetECross(Int_t absID, Double_t tcell, AliVCaloCells *cells, Int_t bc, Float_t cellMinEn=0., Bool_t useWeight=kFALSE, Float_t energyClus=0.)
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
TH2I * GetEMCALChannelStatusMap(Int_t iSM) const
void GetEnergyAndNumberOfCellsInTCard(AliVCluster *clus, Int_t absIdMax, AliVCaloCells *cells, Int_t &nDiff, Int_t &nSame, Float_t &eDiff, Float_t &eSame, Float_t emin=0.)
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:45
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, ESDs.
void SetITSTrackSA(Bool_t isITS)
Bool_t IsExoticCluster(const AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=0)
Int_t GetEMCALL1PhaseInTimeRecalibrationForSM(Int_t iSM, Short_t par=0) const
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.
Short_t fCurrentParNumber
! Current par number
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)
TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const
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)
Short_t fNPars
Number of PARs in run.
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)
Float_t GetEMCALChannelRecalibrationFactor1D(UInt_t iCell) const
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
Bool_t fIsParRun
! flag if run contains PAR
Bool_t GetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow, Int_t &status) const
void RecalculateClusterShowerShapeParameters(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
float Float_t
Definition: External.C:68
Bool_t fUseOuterTrackParam
Use OuterTrackParam not InnerTrackParam, ESDs.
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.
void SetEMCALBadChannelStatusSelection(Bool_t all, Bool_t dead, Bool_t hot, Bool_t warm)
Definition: External.C:228
void RecalibrateCellTimeL1Phase(Int_t iSM, Int_t bc, Double_t &time, Short_t par=0) const
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 GetEMCALChannelStatus1D(Int_t iCell, Int_t &status) const
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.
Bool_t fUse1Drecalib
Flag to use one dimensional recalibration histogram.
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 IsAbsIDsFromTCard(Int_t absId1, Int_t absId2, Int_t &rowDiff, Int_t &colDiff) const
Bool_t IsGoodCluster(AliVCluster *cluster, const AliEMCALGeometry *geom, AliVCaloCells *cells, Int_t bc=-1)
ULong64_t * fGlobalEventID
array with globalID for PAR runs
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)
Bool_t useWeight
Definition: anaTree.C:26
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.
void SetEMCALChannelRecalibrationFactors1D(const TH1S *h)
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
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
void SetEMCALL1PhaseInTimeRecalibrationForSM(Int_t iSM, Int_t c=0, Short_t par=0)
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.
void SetEMCALChannelRecalibrationFactor1D(UInt_t icell, Double_t c=1)
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.
Bool_t fBadStatusSelection[4]
void SetEMCALChannelTimeRecalibrationFactor(Int_t bc, Int_t absID, Double_t c=0, Bool_t isLGon=kFALSE)
Bool_t fCutRequireITSpureSA
ITS pure standalone tracks.
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...
Bool_t fUse1Dmap
Flag to use one 1D bad channel map (for all cells)
void SetEMCALChannelRecalibrationFactors(const TObjArray *map)
Bool_t fUseShaperNonlin
Shaper non linearity correction for towers.
Bool_t fUseTrackDCA
Activate use of aodtrack->GetXYZ or XvYxZv like in AliEMCALRecoUtilsBase::ExtrapolateTrackToEMCalSurf...
Bool_t fLowGain
Switch on or off calibration with low gain channels.
void SetEMCALChannelStatusMap(const TObjArray *map)
Float_t fCutPhi
dPhi cut on matching
Definition: External.C:196
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
Float_t CorrectShaperNonLin(Float_t Emeas, Float_t EcalibHG)
Bool_t fCutEtaPhiSeparate
Cut on dEta and dPhi separately.
Float_t fCutEta
dEta cut on matching
Float_t SmearClusterEnergy(const AliVCluster *clu)
Bool_t fDoUseMergedBC
flag for using one histo for all BCs