AliPhysics  2c6b7ad (2c6b7ad)
AliEMCALRecoUtils.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // ROOT includes
17 #include <TH2F.h>
18 #include <TArrayI.h>
19 #include <TArrayF.h>
20 #include <TObjArray.h>
21 
22 // STEER includes
23 #include "AliVCluster.h"
24 #include "AliVCaloCells.h"
25 #include "AliLog.h"
26 #include "AliPID.h"
27 #include "AliESDEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliESDtrack.h"
30 #include "AliAODTrack.h"
31 #include "AliExternalTrackParam.h"
32 #include "AliESDfriendTrack.h"
33 #include "AliTrackerBase.h"
34 #include "AliMCEvent.h"
35 
36 // EMCAL includes
37 #include "AliEMCALRecoUtils.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliTrackerBase.h"
40 #include "AliEMCALPIDUtils.h"
41 
43 ClassImp(AliEMCALRecoUtils) ;
45 
51 //_____________________________________
53  fParticleType(0), fPosAlgo(0),
54  fW0(0), fShowerShapeCellLocationType(0),
55  fNonLinearityFunction(0), fNonLinearThreshold(0),
56  fSmearClusterEnergy(kFALSE), fRandom(),
57  fCellsRecalibrated(kFALSE), fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
58  fConstantTimeShift(0), fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(), fLowGain(kFALSE),
59  fUseL1PhaseInTimeRecalibration(kFALSE), fEMCALL1PhaseInTimeRecalibration(),
60  fUseRunCorrectionFactors(kFALSE),
61  fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
62  fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
63  fRejectExoticCluster(kFALSE), fRejectExoticCells(kFALSE),
64  fExoticCellFraction(0), fExoticCellDiffTime(0), fExoticCellMinAmplitude(0),
65  fPIDUtils(), fAODFilterMask(0),
66  fAODHybridTracks(0), fAODTPCOnlyTracks(0),
67  fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
68  fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kFALSE), fCutEtaPhiSeparate(kFALSE),
69  fCutR(0), fCutEta(0), fCutPhi(0),
70  fClusterWindow(0), fMass(0),
71  fStepSurface(0), fStepCluster(0),
72  fITSTrackSA(kFALSE), fUseOuterTrackParam(kFALSE), fEMCalSurfaceDistance(440.),
73  fTrackCutsType(0), fCutMinTrackPt(0), fCutMinNClusterTPC(0),
74  fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
75  fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
76  fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0), fCutDCAToVertex2D(kFALSE),
77  fCutRequireITSStandAlone(kFALSE), fCutRequireITSpureSA(kFALSE),
78  fNMCGenerToAccept(0), fMCGenerToAcceptForTrack(1)
79 {
80  // Init parameters
82 
83  for(Int_t j = 0; j < 5; j++) fMCGenerToAccept[j] = "";
84 
85  // Track matching arrays init
88  fResidualPhi = new TArrayF();
89  fResidualEta = new TArrayF();
90  fPIDUtils = new AliEMCALPIDUtils();
91 
92  fBadStatusSelection[0] = kTRUE;
93  fBadStatusSelection[1] = kTRUE;
94  fBadStatusSelection[2] = kTRUE;
95  fBadStatusSelection[3] = kTRUE;
96 }
97 
98 //
99 // Copy constructor.
100 //
101 //______________________________________________________________________
103 : AliEMCALRecoUtilsBase(reco),
112  fLowGain(reco.fLowGain),
117  fEMCALBadChannelMap(NULL),
126  fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0),
127  fResidualPhi( reco.fResidualPhi? new TArrayF(*reco.fResidualPhi):0x0),
129  fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
133  fEMCalSurfaceDistance(440.),
142 {
143  for (Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ;
144  fMisalTransShift[i] = reco.fMisalTransShift[i] ; }
145  for (Int_t i = 0; i < 10 ; i++){ fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
146  for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
147  for (Int_t j = 0; j < 5 ; j++) { fMCGenerToAccept[j] = reco.fMCGenerToAccept[j] ; }
148  for (Int_t j = 0; j < 4 ; j++) { fBadStatusSelection[j] = reco.fBadStatusSelection[j] ; }
149 
150  if(reco.fEMCALBadChannelMap) {
151  // Copy constructor - not taking ownership over calibration histograms
152  fEMCALBadChannelMap = new TObjArray(reco.fEMCALBadChannelMap->GetEntries());
153  fEMCALBadChannelMap->SetOwner(false);
154  for(int ism = 0; ism < reco.fEMCALBadChannelMap->GetEntries(); ism++) fEMCALBadChannelMap->AddAt(reco.fEMCALBadChannelMap->At(ism), ism);
155  }
156 
157  if(reco.fEMCALRecalibrationFactors) {
158  // Copy constructor - not taking ownership over calibration histograms
160  fEMCALRecalibrationFactors->SetOwner(false);
161  for(int ism = 0; ism < reco.fEMCALRecalibrationFactors->GetEntries(); ism++) fEMCALRecalibrationFactors->AddAt(reco.fEMCALRecalibrationFactors->At(ism), ism);
162  }
163 
165  // Copy constructor - not taking ownership over calibration histograms
167  fEMCALTimeRecalibrationFactors->SetOwner(false);
168  for(int ism = 0; ism < reco.fEMCALTimeRecalibrationFactors->GetEntries(); ism++) fEMCALTimeRecalibrationFactors->AddAt(reco.fEMCALTimeRecalibrationFactors->At(ism), ism);
169  }
170 }
171 
175 //______________________________________________________________________
177 {
178  if (this == &reco)return *this;
179  ((TNamed *)this)->operator=(reco);
180 
181  for (Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ;
182  fMisalRotShift[i] = reco.fMisalRotShift[i] ; }
183  for (Int_t i = 0; i < 10 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
184  for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
185 
187  fPosAlgo = reco.fPosAlgo;
188  fW0 = reco.fW0;
190 
194 
197 
200  fLowGain = reco.fLowGain;
201 
204 
206 
209 
212 
218 
219  fPIDUtils = reco.fPIDUtils;
220 
224 
227  fCutR = reco.fCutR;
228  fCutEta = reco.fCutEta;
229  fCutPhi = reco.fCutPhi;
231  fMass = reco.fMass;
232  fStepSurface = reco.fStepSurface;
233  fStepCluster = reco.fStepCluster;
234  fITSTrackSA = reco.fITSTrackSA;
237 
252 
255  for (Int_t j = 0; j < 5 ; j++)
256  fMCGenerToAccept[j] = reco.fMCGenerToAccept[j];
257 
258  //
259  // Assign or copy construct the different TArrays
260  //
261  if (reco.fResidualEta)
262  {
263  if (fResidualEta)
264  *fResidualEta = *reco.fResidualEta;
265  else
266  fResidualEta = new TArrayF(*reco.fResidualEta);
267  }
268  else
269  {
270  if (fResidualEta) delete fResidualEta;
271  fResidualEta = 0;
272  }
273 
274  if (reco.fResidualPhi)
275  {
276  if (fResidualPhi)
277  *fResidualPhi = *reco.fResidualPhi;
278  else
279  fResidualPhi = new TArrayF(*reco.fResidualPhi);
280  }
281  else
282  {
283  if (fResidualPhi) delete fResidualPhi;
284  fResidualPhi = 0;
285  }
286 
287  if (reco.fMatchedTrackIndex)
288  {
289  if (fMatchedTrackIndex)
291  else
293  }
294  else
295  {
297  fMatchedTrackIndex = 0;
298  }
299 
300  if (reco.fMatchedClusterIndex)
301  {
304  else
306  }
307  else
308  {
311  }
312 
313  for (Int_t j = 0; j < 4 ; j++)
315 
317  if(reco.fEMCALBadChannelMap) {
318  // Copy constructor - not taking ownership over calibration histograms
319  fEMCALBadChannelMap = new TObjArray(reco.fEMCALBadChannelMap->GetEntries());
320  fEMCALBadChannelMap->SetOwner(false);
321  for(int ism = 0; ism < reco.fEMCALBadChannelMap->GetEntries(); ism++) fEMCALBadChannelMap->AddAt(reco.fEMCALBadChannelMap->At(ism), ism);
322  }
323 
325  if(reco.fEMCALRecalibrationFactors) {
326  // Copy constructor - not taking ownership over calibration histograms
328  fEMCALRecalibrationFactors->SetOwner(false);
329  for(int ism = 0; ism < reco.fEMCALRecalibrationFactors->GetEntries(); ism++) fEMCALRecalibrationFactors->AddAt(reco.fEMCALRecalibrationFactors->At(ism), ism);
330  }
331 
334  // Copy constructor - not taking ownership over calibration histograms
336  fEMCALTimeRecalibrationFactors->SetOwner(false);
337  for(int ism = 0; ism < reco.fEMCALTimeRecalibrationFactors->GetEntries(); ism++) fEMCALTimeRecalibrationFactors->AddAt(reco.fEMCALTimeRecalibrationFactors->At(ism), ism);
338  }
339 
342  // Copy constructor - not taking ownership over calibration histograms
344  fEMCALL1PhaseInTimeRecalibration->SetOwner(false);
345  for(int ism = 0; ism < reco.fEMCALL1PhaseInTimeRecalibration->GetEntries(); ism++) fEMCALL1PhaseInTimeRecalibration->AddAt(reco.fEMCALL1PhaseInTimeRecalibration->At(ism), ism);
346  }
347 
348  return *this;
349 }
350 
354 //_____________________________________
356 {
358  {
360  }
361 
363  {
365  }
366 
368  {
370  }
371 
372  if (fEMCALBadChannelMap)
373  {
374  delete fEMCALBadChannelMap;
375  }
376 
377  delete fMatchedTrackIndex ;
378  delete fMatchedClusterIndex ;
379  delete fResidualEta ;
380  delete fResidualPhi ;
381  delete fPIDUtils ;
382 
383  InitTrackCuts();
384 }
385 
398 //_______________________________________________________________________________
400  Float_t & amp, Double_t & time,
401  AliVCaloCells* cells)
402 {
403  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
404 
405  if(!geom)
406  {
407  AliError("No instance of the geometry is available");
408  return kFALSE;
409  }
410 
411  if ( absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules() )
412  return kFALSE;
413 
414  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1, status=0;
415 
416  if (!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
417  {
418  // cell absID does not exist
419  amp=0; time = 1.e9;
420  return kFALSE;
421  }
422 
423  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
424 
425  // Do not include bad channels found in analysis,
427  {
428  Bool_t bad = GetEMCALChannelStatus(imod, ieta, iphi,status);
429 
430  if ( status > 0 )
431  AliDebug(1,Form("Channel absId %d, status %d, set as bad %d",absID, status, bad));
432 
433  if ( bad ) return kFALSE;
434  }
435 
436  //Recalibrate energy
437  amp = cells->GetCellAmplitude(absID);
439  amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
440 
441  // Recalibrate time
442  time = cells->GetCellTime(absID);
443  time-=fConstantTimeShift*1e-9; // only in case of old Run1 simulation
444  Bool_t isLowGain = !(cells->GetCellHighGain(absID));//HG = false -> LG = true
445 
446  RecalibrateCellTime(absID,bc,time,isLowGain);
447 
448  //Recalibrate time with L1 phase
449  RecalibrateCellTimeL1Phase(imod, bc, time);
450 
451  return kTRUE;
452 }
453 
463 //_____________________________________________________________________________
465  const AliVCluster* cluster,
466  AliVCaloCells* cells)
467 {
468  if (!cluster)
469  {
470  AliInfo("Cluster pointer null!");
471  return kFALSE;
472  }
473 
474  // If the distance to the border is 0 or negative just exit accept all clusters
475  if (cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 )
476  return kTRUE;
477 
478  Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
479  Bool_t shared = kFALSE;
480  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
481 
482  AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
483  absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
484 
485  if (absIdMax==-1) return kFALSE;
486 
487  // Check if the cell is close to the borders:
488  Bool_t okrow = kFALSE;
489  Bool_t okcol = kFALSE;
490 
491  if (iSM < 0 || iphi < 0 || ieta < 0 )
492  {
493  AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
494  iSM,ieta,iphi));
495  return kFALSE; // trick coverity
496  }
497 
498  // Check rows/phi
499  Int_t iPhiLast = 24;
500  if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_Half ) iPhiLast /= 2;
501  else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_3rd ) iPhiLast /= 3;// 1/3 sm case
502  else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Ext ) iPhiLast /= 3;// 1/3 sm case
503 
504  if(iphi >= fNCellsFromEMCALBorder && iphi < iPhiLast - fNCellsFromEMCALBorder) okrow = kTRUE;
505 
506  // Check columns/eta
507  Int_t iEtaLast = 48;
508  if ( !fNoEMCALBorderAtEta0 || geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Standard )
509  {
510  // consider inner border
511  if ( geom->IsDCALSM(iSM) ) iEtaLast = iEtaLast*2/3;
512 
513  if ( ieta > fNCellsFromEMCALBorder && ieta < iEtaLast-fNCellsFromEMCALBorder ) okcol = kTRUE;
514  }
515  else
516  {
517  if (iSM%2==0)
518  {
519  if (ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
520  }
521  else
522  {
523  if (ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE;
524  }
525  }//eta 0 not checked
526 
527  AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
528  fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
529 
530  if (okcol && okrow)
531  {
532  return kTRUE;
533  }
534  else
535  {
536  AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
537 
538  return kFALSE;
539  }
540 }
541 
552 //_______________________________________________________________________________
554  const UShort_t* cellList,
555  Int_t nCells)
556 {
557  if (!fRemoveBadChannels) return kFALSE;
558  if (!fEMCALBadChannelMap) return kFALSE;
559 
560  Int_t icol = -1;
561  Int_t irow = -1;
562  Int_t imod = -1;
563  for (Int_t iCell = 0; iCell<nCells; iCell++)
564  {
565  //Get the column and row
566  Int_t iTower = -1, iIphi = -1, iIeta = -1;
567  geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
568 
569  if (fEMCALBadChannelMap->GetEntries() <= imod) continue;
570 
571  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
572 
573  Int_t status = 0;
574  if (GetEMCALChannelStatus(imod, icol, irow, status))
575  {
576  AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d, status %d\n",imod, icol, irow, status));
577  return kTRUE;
578  }
579  }// cell cluster loop
580 
581  return kFALSE;
582 }
583 
595 //___________________________________________________________________________
597  AliVCaloCells* cells, Int_t bc)
598 {
599  AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
600 
601  if(!geom)
602  {
603  AliError("No instance of the geometry is available");
604  return -1;
605  }
606 
607  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
608  geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
609  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
610 
611  // Get close cells index, energy and time, not in corners
612 
613  Int_t absID1 = -1;
614  Int_t absID2 = -1;
615 
616  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
617  if ( iphi > 0 ) absID2 = geom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
618 
619  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
620 
621  Int_t absID3 = -1;
622  Int_t absID4 = -1;
623 
624  if ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) )
625  {
626  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
627  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
628  }
629  else if ( ieta == 0 && imod%2 )
630  {
631  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
632  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
633  }
634  else
635  {
636  if ( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
637  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
638  if ( ieta > 0 )
639  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
640  }
641 
642  //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4);
643 
644  Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
645  Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
646 
647  AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
648  AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
649  AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
650  AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
651 
652  if (TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
653  if (TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
654  if (TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
655  if (TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
656 
657  return ecell1+ecell2+ecell3+ecell4;
658 }
659 
670 //_____________________________________________________________________________________________
671 Bool_t AliEMCALRecoUtils::IsExoticCell(Int_t absID, AliVCaloCells* cells, Int_t bc)
672 {
673  if (!fRejectExoticCells) return kFALSE;
674 
675  Float_t ecell = 0;
676  Double_t tcell = 0;
677  Bool_t accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
678 
679  if (!accept) return kTRUE; // reject this cell
680 
681  if (ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
682 
683  Float_t eCross = GetECross(absID,tcell,cells,bc);
684 
685  if (1-eCross/ecell > fExoticCellFraction)
686  {
687  AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",
688  absID,ecell,eCross,1-eCross/ecell));
689  return kTRUE;
690  }
691 
692  return kFALSE;
693 }
694 
704 //___________________________________________________________________
705 Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster,
706  AliVCaloCells *cells,
707  Int_t bc)
708 {
709  if (!cluster)
710  {
711  AliInfo("Cluster pointer null!");
712  return kFALSE;
713  }
714 
715  if (!fRejectExoticCluster) return kFALSE;
716 
717  // Get highest energy tower
718  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
719 
720  if(!geom)
721  {
722  AliError("No instance of the geometry is available");
723  return kFALSE;
724  }
725 
726  Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
727  Bool_t shared = kFALSE;
728  GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
729 
730  return IsExoticCell(absId,cells,bc);
731 }
732 
741 //_______________________________________________________________________
742 Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster)
743 {
744  if (!cluster)
745  {
746  AliInfo("Cluster pointer null!");
747  return 0;
748  }
749 
750  Float_t energy = cluster->E() ;
751  Float_t rdmEnergy = energy ;
752 
753  if (fSmearClusterEnergy)
754  {
755  rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
756  fSmearClusterParam[1] * energy +
757  fSmearClusterParam[2] );
758  AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
759  }
760 
761  return rdmEnergy;
762 }
763 
772 //____________________________________________________________________________
774 {
775  if (!cluster)
776  {
777  AliInfo("Cluster pointer null!");
778  return 0;
779  }
780 
781  Float_t energy = cluster->E();
782  if (energy < 0.100)
783  {
784  // Clusters with less than 100 MeV or negative are not possible
785  AliInfo(Form("Too low cluster energy!, E = %f < 0.100 GeV",energy));
786  return 0;
787  }
788  else if(energy > 300.)
789  {
790  AliInfo(Form("Too high cluster energy!, E = %f GeV, do not apply non linearity",energy));
791  return energy;
792  }
793 
794  switch (fNonLinearityFunction)
795  {
796  case kPi0MC:
797  {
798  //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
799  //fNonLinearityParams[0] = 1.014;
800  //fNonLinearityParams[1] =-0.03329;
801  //fNonLinearityParams[2] =-0.3853;
802  //fNonLinearityParams[3] = 0.5423;
803  //fNonLinearityParams[4] =-0.4335;
804  energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
805  ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
806  exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
807  break;
808  }
809 
810  case kPi0MCv2:
811  {
812  //Non-Linearity correction (from MC with function [0]/((x+[1])^[2]))+1;
813  //fNonLinearityParams[0] = 3.11111e-02;
814  //fNonLinearityParams[1] =-5.71666e-02;
815  //fNonLinearityParams[2] = 5.67995e-01;
816 
817  energy *= fNonLinearityParams[0]/TMath::Power(energy+fNonLinearityParams[1],fNonLinearityParams[2])+1;
818  break;
819  }
820 
821  case kPi0MCv3:
822  {
823  //Same as beam test corrected, change parameters
824  //fNonLinearityParams[0] = 9.81039e-01
825  //fNonLinearityParams[1] = 1.13508e-01;
826  //fNonLinearityParams[2] = 1.00173e+00;
827  //fNonLinearityParams[3] = 9.67998e-02;
828  //fNonLinearityParams[4] = 2.19381e+02;
829  //fNonLinearityParams[5] = 6.31604e+01;
830  //fNonLinearityParams[6] = 1;
831  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
832 
833  break;
834  }
835 
836 
837  case kPi0GammaGamma:
838  {
839  //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
840  //fNonLinearityParams[0] = 1.04;
841  //fNonLinearityParams[1] = -0.1445;
842  //fNonLinearityParams[2] = 1.046;
843  energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
844  break;
845  }
846 
847  case kPi0GammaConversion:
848  {
849  //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
850  //fNonLinearityParams[0] = 0.139393/0.1349766;
851  //fNonLinearityParams[1] = 0.0566186;
852  //fNonLinearityParams[2] = 0.982133;
853  energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
854 
855  break;
856  }
857 
858  case kBeamTest:
859  {
860  //From beam test, Alexei's results, for different ZS thresholds
861  // th=30 MeV; th = 45 MeV; th = 75 MeV
862  //fNonLinearityParams[0] = 1.007; 1.003; 1.002
863  //fNonLinearityParams[1] = 0.894; 0.719; 0.797
864  //fNonLinearityParams[2] = 0.246; 0.334; 0.358
865  //Rescale the param[0] with 1.03
866  energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
867 
868  break;
869  }
870 
871  case kBeamTestCorrected:
872  {
873  // From beam test, corrected for material between beam and EMCAL
874  //fNonLinearityParams[0] = 0.99078
875  //fNonLinearityParams[1] = 0.161499;
876  //fNonLinearityParams[2] = 0.655166;
877  //fNonLinearityParams[3] = 0.134101;
878  //fNonLinearityParams[4] = 163.282;
879  //fNonLinearityParams[5] = 23.6904;
880  //fNonLinearityParams[6] = 0.978;
881  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
882 
883  break;
884  }
885 
887  {
888  // From beam test, corrected for material between beam and EMCAL
889  // Different parametrization to kBeamTestCorrected
890  //fNonLinearityParams[0] = 0.983504;
891  //fNonLinearityParams[1] = 0.210106;
892  //fNonLinearityParams[2] = 0.897274;
893  //fNonLinearityParams[3] = 0.0829064;
894  //fNonLinearityParams[4] = 152.299;
895  //fNonLinearityParams[5] = 31.5028;
896  //fNonLinearityParams[6] = 0.968;
897  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
898 
899  break;
900  }
901 
903  {
904  // Same function as kBeamTestCorrected, different default parametrization.
905  //fNonLinearityParams[0] = 0.976941;
906  //fNonLinearityParams[1] = 0.162310;
907  //fNonLinearityParams[2] = 1.08689;
908  //fNonLinearityParams[3] = 0.0819592;
909  //fNonLinearityParams[4] = 152.338;
910  //fNonLinearityParams[5] = 30.9594;
911  //fNonLinearityParams[6] = 0.9615;
912  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
913 
914  break;
915  }
916 
918  {
919  // New parametrization of kBeamTestCorrected,
920  // fitting new points for E>100 GeV.
921  // I should have same performance as v3 in the low energies
922  // See EMCal meeting 21/09/2018 slides
923  // https://indico.cern.ch/event/759154/contributions/3148448/attachments/1721042/2778585/nonLinearityUpdate.pdf
924  // and jira ticket EMCAL-190
925  // Not very smart copy pasting the same function for each new parametrization, need to think how to do it better.
926 
927 // fNonLinearityParams[0] = 0.9892;
928 // fNonLinearityParams[1] = 0.1976;
929 // fNonLinearityParams[2] = 0.865;
930 // fNonLinearityParams[3] = 0.06775;
931 // fNonLinearityParams[4] = 156.6;
932 // fNonLinearityParams[5] = 47.18;
933 // fNonLinearityParams[6] = 0.97;
934 
935  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
936 
937  break;
938  }
939 
940  case kSDMv5:
941  {
942  //Based on fit to the MC/data using kNoCorrection on the data - utilizes symmetric decay method and kPi0MCv5(MC) - 28 Oct 2013
943  //fNonLinearityParams[0] = 1.0;
944  //fNonLinearityParams[1] = 6.64778e-02;
945  //fNonLinearityParams[2] = 1.570;
946  //fNonLinearityParams[3] = 9.67998e-02;
947  //fNonLinearityParams[4] = 2.19381e+02;
948  //fNonLinearityParams[5] = 6.31604e+01;
949  //fNonLinearityParams[6] = 1.01286;
950  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));
951 
952  break;
953  }
954 
955  case kPi0MCv5:
956  {
957  //Based on comparing MC truth information to the reconstructed energy of clusters.
958  //fNonLinearityParams[0] = 1.0;
959  //fNonLinearityParams[1] = 6.64778e-02;
960  //fNonLinearityParams[2] = 1.570;
961  //fNonLinearityParams[3] = 9.67998e-02;
962  //fNonLinearityParams[4] = 2.19381e+02;
963  //fNonLinearityParams[5] = 6.31604e+01;
964  //fNonLinearityParams[6] = 1.01286;
965  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
966 
967  break;
968  }
969 
970  case kSDMv6:
971  {
972  //Based on fit to the MC/data using kNoCorrection on the data
973  // - utilizes symmetric decay method and kPi0MCv6(MC) - 09 Dec 2014
974  // - parameters constrained by the test beam data as well
975  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
976  //fNonLinearityParams[0] = 1.0;
977  //fNonLinearityParams[1] = 0.237767;
978  //fNonLinearityParams[2] = 0.651203;
979  //fNonLinearityParams[3] = 0.183741;
980  //fNonLinearityParams[4] = 155.427;
981  //fNonLinearityParams[5] = 17.0335;
982  //fNonLinearityParams[6] = 0.987054;
983  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
984 
985  break;
986  }
987 
988  case kPi0MCv6:
989  {
990  //Based on comparing MC truth information to the reconstructed energy of clusters.
991  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
992  //fNonLinearityParams[0] = 1.0;
993  //fNonLinearityParams[1] = 0.0797873;
994  //fNonLinearityParams[2] = 1.68322;
995  //fNonLinearityParams[3] = 0.0806098;
996  //fNonLinearityParams[4] = 244.586;
997  //fNonLinearityParams[5] = 116.938;
998  //fNonLinearityParams[6] = 1.00437;
999  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1000 
1001  break;
1002  }
1003 
1004  case kPCMv1:
1005  {
1006  //based on symmetric decays of pi0 meson
1007  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1008  // parameters vary from MC to MC
1009  //fNonLinearityParams[0] = 0.984876;
1010  //fNonLinearityParams[1] = -9.999609;
1011  //fNonLinearityParams[2] = -4.999891;
1012  //fNonLinearityParams[3] = 0.;
1013  //fNonLinearityParams[4] = 0.;
1014  //fNonLinearityParams[5] = 0.;
1015  //fNonLinearityParams[6] = 0.;
1016  energy /= TMath::Power(fNonLinearityParams[0] + exp(fNonLinearityParams[1] + fNonLinearityParams[2]*energy),2);//result coming from calo-conv method
1017 
1018  break;
1019  }
1020 
1021  case kPCMplusBTCv1:
1022  {
1023  //convolution of TestBeamCorrectedv3 with PCM method
1024  //Based on comparing MC truth information to the reconstructed energy of clusters.
1025  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1026  // parameters vary from MC to MC
1027  //fNonLinearityParams[0] = 0.976941;
1028  //fNonLinearityParams[1] = 0.162310;
1029  //fNonLinearityParams[2] = 1.08689;
1030  //fNonLinearityParams[3] = 0.0819592;
1031  //fNonLinearityParams[4] = 152.338;
1032  //fNonLinearityParams[5] = 30.9594;
1033  //fNonLinearityParams[6] = 0.9615;
1034  //fNonLinearityParams[7] = 0.984876;
1035  //fNonLinearityParams[8] = -9.999609;
1036  //fNonLinearityParams[9] = -4.999891;
1037  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1038  energy /= TMath::Power(fNonLinearityParams[7] + exp(fNonLinearityParams[8] + fNonLinearityParams[9]*energy),2);//result coming from calo-conv method
1039 
1040  break;
1041  }
1042 
1043  case kPCMsysv1:
1044  {
1045  // Systematic variation of kPCMv1
1046  //Based on comparing MC truth information to the reconstructed energy of clusters.
1047  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1048  // parameters vary from MC to MC
1049  //fNonLinearityParams[0] = 0.0;
1050  //fNonLinearityParams[1] = 1.0;
1051  //fNonLinearityParams[2] = 1.0;
1052  //fNonLinearityParams[3] = 0.0;
1053  //fNonLinearityParams[4] = 1.0;
1054  //fNonLinearityParams[5] = 0.0;
1055  //fNonLinearityParams[6] = 0.0;
1056  energy /= TMath::Power( (fNonLinearityParams[0] + fNonLinearityParams[1] * TMath::Power(energy,fNonLinearityParams[2]) ) /
1057  (fNonLinearityParams[3] + fNonLinearityParams[4] * TMath::Power(energy,fNonLinearityParams[5]) ) + fNonLinearityParams[6], 2);//result coming from calo-conv method
1058 
1059  break;
1060  }
1061 
1062 
1063  case kNoCorrection:
1064  AliDebug(2,"No correction on the energy\n");
1065  break;
1066 
1067  }
1068 
1069  return energy;
1070 }
1071 
1076 //__________________________________________________
1078 {
1079  if (fNonLinearityFunction == kPi0MC) {
1080  fNonLinearityParams[0] = 1.014;
1081  fNonLinearityParams[1] = -0.03329;
1082  fNonLinearityParams[2] = -0.3853;
1083  fNonLinearityParams[3] = 0.5423;
1084  fNonLinearityParams[4] = -0.4335;
1085  }
1086 
1088  fNonLinearityParams[0] = 3.11111e-02;
1089  fNonLinearityParams[1] =-5.71666e-02;
1090  fNonLinearityParams[2] = 5.67995e-01;
1091  }
1092 
1094  fNonLinearityParams[0] = 9.81039e-01;
1095  fNonLinearityParams[1] = 1.13508e-01;
1096  fNonLinearityParams[2] = 1.00173e+00;
1097  fNonLinearityParams[3] = 9.67998e-02;
1098  fNonLinearityParams[4] = 2.19381e+02;
1099  fNonLinearityParams[5] = 6.31604e+01;
1100  fNonLinearityParams[6] = 1;
1101  }
1102 
1104  fNonLinearityParams[0] = 1.04;
1105  fNonLinearityParams[1] = -0.1445;
1106  fNonLinearityParams[2] = 1.046;
1107  }
1108 
1110  fNonLinearityParams[0] = 0.139393;
1111  fNonLinearityParams[1] = 0.0566186;
1112  fNonLinearityParams[2] = 0.982133;
1113  }
1114 
1116  if (fNonLinearThreshold == 30) {
1117  fNonLinearityParams[0] = 1.007;
1118  fNonLinearityParams[1] = 0.894;
1119  fNonLinearityParams[2] = 0.246;
1120  }
1121  if (fNonLinearThreshold == 45) {
1122  fNonLinearityParams[0] = 1.003;
1123  fNonLinearityParams[1] = 0.719;
1124  fNonLinearityParams[2] = 0.334;
1125  }
1126  if (fNonLinearThreshold == 75) {
1127  fNonLinearityParams[0] = 1.002;
1128  fNonLinearityParams[1] = 0.797;
1129  fNonLinearityParams[2] = 0.358;
1130  }
1131  }
1132 
1134  fNonLinearityParams[0] = 0.99078;
1135  fNonLinearityParams[1] = 0.161499;
1136  fNonLinearityParams[2] = 0.655166;
1137  fNonLinearityParams[3] = 0.134101;
1138  fNonLinearityParams[4] = 163.282;
1139  fNonLinearityParams[5] = 23.6904;
1140  fNonLinearityParams[6] = 0.978;
1141  }
1142 
1144  // Parameters until November 2015, use now kBeamTestCorrectedv3
1145  fNonLinearityParams[0] = 0.983504;
1146  fNonLinearityParams[1] = 0.210106;
1147  fNonLinearityParams[2] = 0.897274;
1148  fNonLinearityParams[3] = 0.0829064;
1149  fNonLinearityParams[4] = 152.299;
1150  fNonLinearityParams[5] = 31.5028;
1151  fNonLinearityParams[6] = 0.968;
1152  }
1153 
1155 
1156  // New parametrization of kBeamTestCorrected
1157  // excluding point at 0.5 GeV from Beam Test Data
1158  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1159 
1160  fNonLinearityParams[0] = 0.976941;
1161  fNonLinearityParams[1] = 0.162310;
1162  fNonLinearityParams[2] = 1.08689;
1163  fNonLinearityParams[3] = 0.0819592;
1164  fNonLinearityParams[4] = 152.338;
1165  fNonLinearityParams[5] = 30.9594;
1166  fNonLinearityParams[6] = 0.9615;
1167  }
1168 
1170 
1171  // New parametrization of kBeamTestCorrected,
1172  // fitting new points for E>100 GeV.
1173  // I should have same performance as v3 in the low energies
1174  // See EMCal meeting 21/09/2018 slides
1175  // https://indico.cern.ch/event/759154/contributions/3148448/attachments/1721042/2778585/nonLinearityUpdate.pdf
1176  // and jira ticket EMCAL-190
1177 
1178  fNonLinearityParams[0] = 0.9892;
1179  fNonLinearityParams[1] = 0.1976;
1180  fNonLinearityParams[2] = 0.865;
1181  fNonLinearityParams[3] = 0.06775;
1182  fNonLinearityParams[4] = 156.6;
1183  fNonLinearityParams[5] = 47.18;
1184  fNonLinearityParams[6] = 0.97;
1185  }
1186 
1187  if (fNonLinearityFunction == kSDMv5) {
1188  fNonLinearityParams[0] = 1.0;
1189  fNonLinearityParams[1] = 6.64778e-02;
1190  fNonLinearityParams[2] = 1.570;
1191  fNonLinearityParams[3] = 9.67998e-02;
1192  fNonLinearityParams[4] = 2.19381e+02;
1193  fNonLinearityParams[5] = 6.31604e+01;
1194  fNonLinearityParams[6] = 1.01286;
1195  }
1196 
1198  fNonLinearityParams[0] = 1.0;
1199  fNonLinearityParams[1] = 6.64778e-02;
1200  fNonLinearityParams[2] = 1.570;
1201  fNonLinearityParams[3] = 9.67998e-02;
1202  fNonLinearityParams[4] = 2.19381e+02;
1203  fNonLinearityParams[5] = 6.31604e+01;
1204  fNonLinearityParams[6] = 1.01286;
1205  }
1206 
1207  if (fNonLinearityFunction == kSDMv6) {
1208  fNonLinearityParams[0] = 1.0;
1209  fNonLinearityParams[1] = 0.237767;
1210  fNonLinearityParams[2] = 0.651203;
1211  fNonLinearityParams[3] = 0.183741;
1212  fNonLinearityParams[4] = 155.427;
1213  fNonLinearityParams[5] = 17.0335;
1214  fNonLinearityParams[6] = 0.987054;
1215  }
1216 
1218  fNonLinearityParams[0] = 1.0;
1219  fNonLinearityParams[1] = 0.0797873;
1220  fNonLinearityParams[2] = 1.68322;
1221  fNonLinearityParams[3] = 0.0806098;
1222  fNonLinearityParams[4] = 244.586;
1223  fNonLinearityParams[5] = 116.938;
1224  fNonLinearityParams[6] = 1.00437;
1225  }
1226 
1227 if (fNonLinearityFunction == kPCMv1) {
1228  //parameters change from MC production to MC production, they need to set for each period
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  }
1237 
1239  // test beam corrected values convoluted with symmetric meson decays values
1240  // for test beam:
1241  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1242  // for PCM method:
1243  // https://aliceinfo.cern.ch/Notes/node/211
1244  fNonLinearityParams[0] = 0.976941;
1245  fNonLinearityParams[1] = 0.162310;
1246  fNonLinearityParams[2] = 1.08689;
1247  fNonLinearityParams[3] = 0.0819592;
1248  fNonLinearityParams[4] = 152.338;
1249  fNonLinearityParams[5] = 30.9594;
1250  fNonLinearityParams[6] = 0.9615;
1251  fNonLinearityParams[7] = 0.984876;
1252  fNonLinearityParams[8] = -9.999609;
1253  fNonLinearityParams[9] = -4.999891;
1254  }
1255 
1257  //systematics for kPCMv1
1258  // for PCM method:
1259  // https://aliceinfo.cern.ch/Notes/node/211
1260  fNonLinearityParams[0] = 0.0;
1261  fNonLinearityParams[1] = 1.0;
1262  fNonLinearityParams[2] = 1.0;
1263  fNonLinearityParams[3] = 0.0;
1264  fNonLinearityParams[4] = 1.0;
1265  fNonLinearityParams[5] = 0.0;
1266  fNonLinearityParams[6] = 0.0;
1267  }
1268 }
1269 
1277 //____________________________________________________________________
1279 {
1280  fBadStatusSelection[0] = all; // declare all as bad if true, never mind the other settings
1281  fBadStatusSelection[1] = dead;
1282  fBadStatusSelection[2] = hot;
1283  fBadStatusSelection[3] = warm;
1284 }
1285 
1296 //____________________________________________________________________
1298 {
1299  if(fEMCALBadChannelMap)
1300  status = (Int_t) ((TH2I*)fEMCALBadChannelMap->At(iSM))->GetBinContent(iCol,iRow);
1301  else
1302  status = 0; // Channel is ok by default
1303 
1304  if ( status == AliCaloCalibPedestal::kAlive )
1305  {
1306  return kFALSE; // Good channel
1307  }
1308  else
1309  {
1310  if ( fBadStatusSelection[0] == kTRUE )
1311  {
1312  return kTRUE; // consider bad hot, dead and warm
1313  }
1314  else
1315  {
1316  if ( fBadStatusSelection[AliCaloCalibPedestal::kDead] == kTRUE &&
1317  status == AliCaloCalibPedestal::kDead )
1318  return kTRUE; // consider bad dead
1319  else if ( fBadStatusSelection[AliCaloCalibPedestal::kHot] == kTRUE &&
1320  status == AliCaloCalibPedestal::kHot )
1321  return kTRUE; // consider bad hot
1322  else if ( fBadStatusSelection[AliCaloCalibPedestal::kWarning] == kTRUE &&
1323  status == AliCaloCalibPedestal::kWarning )
1324  return kTRUE; // consider bad warm
1325  }
1326  }
1327 
1328  AliWarning(Form("Careful, bad channel selection not properly done: ism %d, icol %d, irow %d, status %d,\n"
1329  " fBadAll %d, fBadHot %d, fBadWarm %d, fBadDead %d",
1330  iSM, iCol, iRow, status,
1333 
1334  return kFALSE; // if everything fails, accept it.
1335 }
1336 
1350 //____________________________________________________________________
1351 void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
1352  AliVCaloCells* cells,
1353  const AliVCluster* clu,
1354  Int_t & absId,
1355  Int_t & iSupMod,
1356  Int_t & ieta,
1357  Int_t & iphi,
1358  Bool_t & shared)
1359 {
1360  Double_t eMax = -1.;
1361  Double_t eCell = -1.;
1362  Float_t fraction = 1.;
1363  Float_t recalFactor = 1.;
1364  Int_t cellAbsId = -1 ;
1365 
1366  Int_t iTower = -1;
1367  Int_t iIphi = -1;
1368  Int_t iIeta = -1;
1369  Int_t iSupMod0= -1;
1370 
1371  if (!clu)
1372  {
1373  AliInfo("Cluster pointer null!");
1374  absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
1375  return;
1376  }
1377 
1378  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1379  {
1380  cellAbsId = clu->GetCellAbsId(iDig);
1381  fraction = clu->GetCellAmplitudeFraction(iDig);
1382  //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
1383  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1384 
1385  geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1386  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1387 
1388  if (iDig==0)
1389  {
1390  iSupMod0=iSupMod;
1391  } else if (iSupMod0!=iSupMod)
1392  {
1393  shared = kTRUE;
1394  //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
1395  }
1396 
1398  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1399 
1400  eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
1401  //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
1402  if (eCell > eMax)
1403  {
1404  eMax = eCell;
1405  absId = cellAbsId;
1406  //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
1407  }
1408  }// cell loop
1409 
1410  //Get from the absid the supermodule, tower and eta/phi numbers
1411  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1412  //Gives SuperModule and Tower numbers
1413  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1414  //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
1415  //printf("Max end---\n");
1416 }
1417 
1426 //_________________________________________________________
1428 {
1429  if (eCell > 0 && eCluster > 0)
1430  {
1431  if ( fW0 > 0 ) return TMath::Max( 0., fW0 + TMath::Log( eCell / eCluster ) ) ;
1432  else return TMath::Log( eCluster / eCell ) ;
1433  }
1434  else
1435  return 0. ;
1436 }
1437 
1441 //______________________________________
1443 {
1444  fParticleType = kPhoton;
1445  fPosAlgo = kUnchanged;
1446  fW0 = 4.5;
1447 
1449  fNonLinearThreshold = 30;
1450 
1451  fExoticCellFraction = 0.97;
1452  fExoticCellDiffTime = 1e6;
1454 
1455  fAODFilterMask = 128;
1456  fAODHybridTracks = kFALSE;
1457  fAODTPCOnlyTracks = kTRUE;
1458 
1459  fCutEtaPhiSum = kTRUE;
1460  fCutEtaPhiSeparate = kFALSE;
1461 
1462  fCutR = 0.05;
1463  fCutEta = 0.025;
1464  fCutPhi = 0.05;
1465 
1466  fClusterWindow = 100;
1467  fMass = 0.139;
1468 
1469  fStepSurface = 20.;
1470  fStepCluster = 5.;
1472 
1473  fCutMinTrackPt = 0;
1474  fCutMinNClusterTPC = -1;
1475  fCutMinNClusterITS = -1;
1476 
1477  fCutMaxChi2PerClusterTPC = 1e10;
1478  fCutMaxChi2PerClusterITS = 1e10;
1479 
1480  fCutRequireTPCRefit = kFALSE;
1481  fCutRequireITSRefit = kFALSE;
1482  fCutAcceptKinkDaughters = kFALSE;
1483 
1484  fCutMaxDCAToVertexXY = 1e10;
1485  fCutMaxDCAToVertexZ = 1e10;
1486  fCutDCAToVertex2D = kFALSE;
1487 
1488  fCutRequireITSStandAlone = kFALSE; //MARCEL
1489  fCutRequireITSpureSA = kFALSE; //Marcel
1490 
1491  // Misalignment matrices
1492  for (Int_t i = 0; i < 15 ; i++)
1493  {
1494  fMisalTransShift[i] = 0.;
1495  fMisalRotShift[i] = 0.;
1496  }
1497 
1498  // Non linearity
1499  for (Int_t i = 0; i < 10 ; i++) fNonLinearityParams[i] = 0.;
1500 
1501  // For kBeamTestCorrectedv2 case, but default is no correction
1502  fNonLinearityParams[0] = 0.983504;
1503  fNonLinearityParams[1] = 0.210106;
1504  fNonLinearityParams[2] = 0.897274;
1505  fNonLinearityParams[3] = 0.0829064;
1506  fNonLinearityParams[4] = 152.299;
1507  fNonLinearityParams[5] = 31.5028;
1508  fNonLinearityParams[6] = 0.968;
1509 
1510  // Cluster energy smearing
1511  fSmearClusterEnergy = kFALSE;
1512  fSmearClusterParam[0] = 0.07; // * sqrt E term
1513  fSmearClusterParam[1] = 0.00; // * E term
1514  fSmearClusterParam[2] = 0.00; // constant
1515 }
1516 
1520 //_____________________________________________________
1522 {
1523  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1524 
1525  // In order to avoid rewriting the same histograms
1526  Bool_t oldStatus = TH1::AddDirectoryStatus();
1527  TH1::AddDirectory(kFALSE);
1528 
1530  for (int i = 0; i < 22; i++)
1531  fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1532  Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
1533  //Init the histograms with 1
1534  for (Int_t sm = 0; sm < 22; sm++)
1535  {
1536  for (Int_t i = 0; i < 48; i++)
1537  {
1538  for (Int_t j = 0; j < 24; j++)
1539  {
1541  }
1542  }
1543  }
1544 
1545  fEMCALRecalibrationFactors->SetOwner(kTRUE);
1546  fEMCALRecalibrationFactors->Compress();
1547 
1548  // In order to avoid rewriting the same histograms
1549  TH1::AddDirectory(oldStatus);
1550 }
1551 
1555 //_________________________________________________________
1557 {
1558  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1559 
1560  // In order to avoid rewriting the same histograms
1561  Bool_t oldStatus = TH1::AddDirectoryStatus();
1562  TH1::AddDirectory(kFALSE);
1563 
1566 
1567  for (int i = 0; i < 4; i++)
1568  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1569  Form("hAllTimeAvBC%d",i),
1570  48*24*22,0.,48*24*22) );
1571  // Init the histograms with 0
1572  for (Int_t iBC = 0; iBC < 4; iBC++)
1573  {
1574  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1575  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kFALSE);
1576  }
1577 
1578  if(fLowGain) {
1579  for (int iBC = 0; iBC < 4; iBC++) {
1580  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvLGBC%d",iBC),
1581  Form("hAllTimeAvLGBC%d",iBC),
1582  48*24*22,0.,48*24*22) );
1583  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1584  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kTRUE);
1585  }
1586  }
1587 
1588 
1589 
1590  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1591  fEMCALTimeRecalibrationFactors->Compress();
1592 
1593  // In order to avoid rewriting the same histograms
1594  TH1::AddDirectory(oldStatus);
1595 }
1596 
1600 //____________________________________________________
1602 {
1603  AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1604 
1605  // In order to avoid rewriting the same histograms
1606  Bool_t oldStatus = TH1::AddDirectoryStatus();
1607  TH1::AddDirectory(kFALSE);
1608 
1609  fEMCALBadChannelMap = new TObjArray(22);
1610  //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
1611 
1612  for (int i = 0; i < 22; i++)
1613  fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1614 
1615  fEMCALBadChannelMap->SetOwner(kTRUE);
1616  fEMCALBadChannelMap->Compress();
1617 
1618  // In order to avoid rewriting the same histograms
1619  TH1::AddDirectory(oldStatus);
1620 }
1621 
1625 //___________________________________________________________
1627 {
1628  AliDebug(2,"AliEMCALRecoUtils::InitEMCALL1PhaseInTimeRecalibrationFactors()");
1629 
1630  // In order to avoid rewriting the same histograms
1631  Bool_t oldStatus = TH1::AddDirectoryStatus();
1632  TH1::AddDirectory(kFALSE);
1633 
1635 
1636  fEMCALL1PhaseInTimeRecalibration->Add(new TH1C("h0","EMCALL1phaseForSM", 22, 0, 22));
1637 
1638  for (Int_t i = 0; i < 22; i++) //loop over SMs, default value = 0
1640 
1641  fEMCALL1PhaseInTimeRecalibration->SetOwner(kTRUE);
1643 
1644  // In order to avoid rewriting the same histograms
1645  TH1::AddDirectory(oldStatus);
1646 }
1647 
1657 //____________________________________________________________________________
1658 void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
1659  AliVCluster * cluster,
1660  AliVCaloCells * cells,
1661  Int_t bc)
1662 {
1663  if (!cluster)
1664  {
1665  AliInfo("Cluster pointer null!");
1666  return;
1667  }
1668 
1669  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1670  UShort_t * index = cluster->GetCellsAbsId() ;
1671  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1672  Int_t ncells = cluster->GetNCells();
1673 
1674  // Initialize some used variables
1675  Float_t energy = 0;
1676  Int_t absId =-1;
1677  Int_t icol =-1, irow =-1, imod=1;
1678  Float_t factor = 1, frac = 0;
1679  Int_t absIdMax = -1;
1680  Float_t emax = 0;
1681 
1682  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1683  for (Int_t icell = 0; icell < ncells; icell++)
1684  {
1685  absId = index[icell];
1686  frac = fraction[icell];
1687  if (frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
1688 
1690  {
1691  // Energy
1692  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1693  geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1694  if (fEMCALRecalibrationFactors->GetEntries() <= imod)
1695  continue;
1696  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1697  factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1698 
1699  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1700  imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1701 
1702  }
1703 
1704  energy += cells->GetCellAmplitude(absId)*factor*frac;
1705 
1706  if (emax < cells->GetCellAmplitude(absId)*factor*frac)
1707  {
1708  emax = cells->GetCellAmplitude(absId)*factor*frac;
1709  absIdMax = absId;
1710  }
1711  }
1712 
1713  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
1714 
1715  cluster->SetE(energy);
1716 
1717  // Recalculate time of cluster
1718  Double_t timeorg = cluster->GetTOF();
1719  Bool_t isLowGain = !(cells->GetCellHighGain(absIdMax));//HG = false -> LG = true
1720 
1721  Double_t time = cells->GetCellTime(absIdMax);
1723  RecalibrateCellTime(absIdMax,bc,time,isLowGain);
1724  time-=fConstantTimeShift*1e-9; // only in case of Run1 old simulations
1725 
1726  // Recalibrate time with L1 phase
1728  RecalibrateCellTimeL1Phase(imod, bc, time);
1729 
1730  cluster->SetTOF(time);
1731 
1732  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
1733 }
1734 
1742 //_______________________________________________________________________
1743 void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, Int_t bc)
1744 {
1746  return;
1747 
1748  if (!cells)
1749  {
1750  AliInfo("Cells pointer null!");
1751  return;
1752  }
1753 
1754  Short_t absId =-1;
1755  Bool_t accept = kFALSE;
1756  Float_t ecell = 0;
1757  Double_t tcell = 0;
1758  Double_t ecellin = 0;
1759  Double_t tcellin = 0;
1760  Int_t mclabel = -1;
1761  Double_t efrac = 0;
1762 
1763  Int_t nEMcell = cells->GetNumberOfCells() ;
1764  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1765  {
1766  cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
1767 
1768  accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
1769  if (!accept)
1770  {
1771  ecell = 0;
1772  tcell = -1;
1773  }
1774 
1775  // Set new values
1776  cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
1777  }
1778 
1779  fCellsRecalibrated = kTRUE;
1780 }
1781 
1790 //_______________________________________________________________________________________________________
1791 void AliEMCALRecoUtils::RecalibrateCellTime(Int_t absId, Int_t bc, Double_t & celltime, Bool_t isLGon) const
1792 {
1793  if (!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0) {
1794  if(fLowGain)
1795  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,isLGon)*1.e-9;
1796  else
1797  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,kFALSE)*1.e-9;
1798  }
1799 }
1800 
1808 //_______________________________________________________________________________________________________
1810 {
1811  if (!fCellsRecalibrated && IsL1PhaseInTimeRecalibrationOn() && bc >= 0)
1812  {
1813  bc=bc%4;
1814 
1815  Float_t offsetPerSM=0.;
1816  Int_t l1PhaseShift = GetEMCALL1PhaseInTimeRecalibrationForSM(iSM);
1817  Int_t l1Phase=l1PhaseShift & 3; //bit operation
1818 
1819  if(bc >= l1Phase)
1820  offsetPerSM = (bc - l1Phase)*25;
1821  else
1822  offsetPerSM = (bc - l1Phase + 4)*25;
1823 
1824  Int_t l1shiftOffset=l1PhaseShift>>2; //bit operation
1825  l1shiftOffset*=25;
1826 
1827  celltime -= offsetPerSM*1.e-9;
1828  celltime -= l1shiftOffset*1.e-9;
1829  }
1830 }
1831 
1847 //______________________________________________________________________________
1848 void AliEMCALRecoUtils::RecalculateCellLabelsRemoveAddedGenerator( Int_t absID, AliVCluster* clus, AliMCEvent* mc,
1849  Float_t & amp, TArrayI & labeArr, TArrayF & eDepArr) const
1850 {
1851  TString genName ;
1852  Float_t eDepFrac[4];
1853 
1854  Float_t edepTotFrac = 1;
1855  Bool_t found = kFALSE;
1856  Float_t ampOrg = amp;
1857 
1858  //
1859  // Get the energy deposition fraction from cluster.
1860  //
1861  for(Int_t icluscell = 0; icluscell < clus->GetNCells(); icluscell++ )
1862  {
1863  if(absID == clus->GetCellAbsId(icluscell))
1864  {
1865  clus->GetCellMCEdepFractionArray(icluscell,eDepFrac);
1866 
1867  found = kTRUE;
1868 
1869  break;
1870  }
1871  }
1872 
1873  if ( !found )
1874  {
1875  AliWarning(Form("Cell abs ID %d NOT found in cluster",absID));
1876  return;
1877  }
1878 
1879  //
1880  // Check if there is a particle contribution from a given generator name.
1881  // If it is not one of the selected generators,
1882  // remove the constribution from the cell.
1883  //
1884  if ( mc && fNMCGenerToAccept > 0 )
1885  {
1886  //printf("Accept contribution from generator? \n");
1887  for(UInt_t imc = 0; imc < 4; imc++)
1888  {
1889  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc )
1890  {
1891  mc->GetCocktailGenerator(clus->GetLabelAt(imc),genName);
1892 
1893  Bool_t generOK = kFALSE;
1894  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
1895  {
1896  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
1897  }
1898 
1899  if ( !generOK )
1900  {
1901  amp-=ampOrg*eDepFrac[imc];
1902 
1903  edepTotFrac-=eDepFrac[imc];
1904 
1905  eDepFrac[imc] = 0;
1906  }
1907 
1908  } // eDep > 0
1909  } // 4 possible loop
1910 
1911  } // accept at least one kind of generator
1912 
1913  //
1914  // Add MC label and Edep to corresponding array (to be used later in digits)
1915  //
1916  Int_t nLabels = 0;
1917  for(UInt_t imc = 0; imc < 4; imc++)
1918  {
1919  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc && edepTotFrac > 0 )
1920  {
1921  nLabels++;
1922 
1923  labeArr.Set(nLabels);
1924  labeArr.AddAt(clus->GetLabelAt(imc), nLabels-1);
1925 
1926  eDepArr.Set(nLabels);
1927  eDepArr.AddAt( (eDepFrac[imc]/edepTotFrac) * amp, nLabels-1);
1928  // use as deposited energy a fraction of the simulated energy (smeared and with noise)
1929  } // edep > 0
1930 
1931  } // mc cell label loop
1932 
1933  //
1934  // If no label found, reject cell
1935  // It can happen to have this case (4 MC labels per cell is not enough for some cases)
1936  // better to remove. To be treated carefully.
1937  //
1938  if ( nLabels == 0 ) amp = 0;
1939 
1940 }
1941 
1949 //______________________________________________________________________________
1950 void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
1951  AliVCaloCells* cells,
1952  AliVCluster* clu)
1953 {
1954  if (!clu)
1955  {
1956  AliInfo("Cluster pointer null!");
1957  return;
1958  }
1959 
1961  else if (fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1962  else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
1963 }
1964 
1973 //_____________________________________________________________________________________________
1975  AliVCaloCells* cells,
1976  AliVCluster* clu)
1977 {
1978  Double_t eCell = 0.;
1979  Float_t fraction = 1.;
1980  Float_t recalFactor = 1.;
1981 
1982  Int_t absId = -1;
1983  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1984  Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
1985  Float_t weight = 0., totalWeight=0.;
1986  Float_t newPos[3] = {-1.,-1.,-1.};
1987  Double_t pLocal[3], pGlobal[3];
1988  Bool_t shared = kFALSE;
1989 
1990  Float_t clEnergy = clu->E(); //Energy already recalibrated previously
1991  if (clEnergy <= 0) return;
1992 
1993  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
1994  Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1995 
1996  //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1997 
1998  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1999  {
2000  absId = clu->GetCellAbsId(iDig);
2001  fraction = clu->GetCellAmplitudeFraction(iDig);
2002  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2003 
2004  if (!fCellsRecalibrated)
2005  {
2006  geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
2007  geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
2008  if (IsRecalibrationOn()) {
2009  recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
2010  }
2011  }
2012 
2013  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2014 
2015  weight = GetCellWeight(eCell,clEnergy);
2016  totalWeight += weight;
2017 
2018  geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
2019  //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
2020  geom->GetGlobal(pLocal,pGlobal,iSupModMax);
2021  //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
2022 
2023  for (int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
2024  }// cell loop
2025 
2026  if (totalWeight>0)
2027  {
2028  for (int i=0; i<3; i++ ) newPos[i] /= totalWeight;
2029  }
2030 
2031  //Float_t pos[]={0,0,0};
2032  //clu->GetPosition(pos);
2033  //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
2034  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
2035 
2036  if (iSupModMax > 1) { //sector 1
2037  newPos[0] +=fMisalTransShift[3];//-=3.093;
2038  newPos[1] +=fMisalTransShift[4];//+=6.82;
2039  newPos[2] +=fMisalTransShift[5];//+=1.635;
2040  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
2041  } else { //sector 0
2042  newPos[0] +=fMisalTransShift[0];//+=1.134;
2043  newPos[1] +=fMisalTransShift[1];//+=8.2;
2044  newPos[2] +=fMisalTransShift[2];//+=1.197;
2045  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
2046  }
2047  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
2048 
2049  clu->SetPosition(newPos);
2050 }
2051 
2060 //____________________________________________________________________________________________
2062  AliVCaloCells* cells,
2063  AliVCluster* clu)
2064 {
2065  Double_t eCell = 1.;
2066  Float_t fraction = 1.;
2067  Float_t recalFactor = 1.;
2068 
2069  Int_t absId = -1;
2070  Int_t iTower = -1;
2071  Int_t iIphi = -1, iIeta = -1;
2072  Int_t iSupMod = -1, iSupModMax = -1;
2073  Int_t iphi = -1, ieta =-1;
2074  Bool_t shared = kFALSE;
2075 
2076  Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
2077 
2078  if (clEnergy <= 0)
2079  return;
2080  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
2081  Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
2082 
2083  Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
2084  Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
2085  Int_t startingSM = -1;
2086 
2087  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
2088  {
2089  absId = clu->GetCellAbsId(iDig);
2090  fraction = clu->GetCellAmplitudeFraction(iDig);
2091  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2092 
2093  if (iDig==0) startingSM = iSupMod;
2094  else if (iSupMod != startingSM) areInSameSM = kFALSE;
2095 
2096  eCell = cells->GetCellAmplitude(absId);
2097 
2098  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2099  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
2100 
2101  if (!fCellsRecalibrated)
2102  {
2103  if (IsRecalibrationOn()) {
2104  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2105  }
2106  }
2107 
2108  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2109 
2110  weight = GetCellWeight(eCell,clEnergy);
2111  if (weight < 0) weight = 0;
2112  totalWeight += weight;
2113  weightedCol += ieta*weight;
2114  weightedRow += iphi*weight;
2115 
2116  //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
2117  }// cell loop
2118 
2119  Float_t xyzNew[]={-1.,-1.,-1.};
2120  if (areInSameSM == kTRUE)
2121  {
2122  //printf("In Same SM\n");
2123  weightedCol = weightedCol/totalWeight;
2124  weightedRow = weightedRow/totalWeight;
2125  geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2126  }
2127  else
2128  {
2129  //printf("In Different SM\n");
2130  geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2131  }
2132 
2133  clu->SetPosition(xyzNew);
2134 }
2135 
2143 //___________________________________________________________________________________________
2145  AliVCaloCells* cells,
2146  AliVCluster * cluster)
2147 {
2148  if (!fRecalDistToBadChannels) return;
2149 
2150  if (!cluster)
2151  {
2152  AliInfo("Cluster pointer null!");
2153  return;
2154  }
2155 
2156  // Get channels map of the supermodule where the cluster is.
2157  Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
2158  Bool_t shared = kFALSE;
2159  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
2160  TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
2161 
2162  Int_t dRrow, dRcol;
2163  Float_t minDist = 10000.;
2164  Float_t dist = 0.;
2165 
2166  // Loop on tower status map
2167  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2168  {
2169  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2170  {
2171  // Check if tower is bad.
2172  if (hMap->GetBinContent(icol,irow)==0) continue;
2173  //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
2174  // iSupMod,icol, irow, icolM,irowM);
2175 
2176  dRrow=TMath::Abs(irowM-irow);
2177  dRcol=TMath::Abs(icolM-icol);
2178  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2179  if (dist < minDist)
2180  {
2181  //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
2182  minDist = dist;
2183  }
2184  }
2185  }
2186 
2187  // In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
2188  if (shared)
2189  {
2190  TH2D* hMap2 = 0;
2191  Int_t iSupMod2 = -1;
2192 
2193  // The only possible combinations are (0,1), (2,3) ... (8,9)
2194  if (iSupMod%2) iSupMod2 = iSupMod-1;
2195  else iSupMod2 = iSupMod+1;
2196  hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
2197 
2198  // Loop on tower status map of second super module
2199  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2200  {
2201  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2202  {
2203  // Check if tower is bad.
2204  if (hMap2->GetBinContent(icol,irow)==0) continue;
2205 
2206  //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",
2207  // iSupMod2,icol, irow,iSupMod,icolM,irowM);
2208  dRrow=TMath::Abs(irow-irowM);
2209 
2210  if (iSupMod%2)
2211  dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
2212  else
2213  dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
2214 
2215  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2216  if (dist < minDist) minDist = dist;
2217  }
2218  }
2219  }// shared cluster in 2 SuperModules
2220 
2221  AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
2222  cluster->SetDistanceToBadChannel(minDist);
2223 }
2224 
2231 //__________________________________________________________________
2232 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
2233 {
2234  if (!cluster)
2235  {
2236  AliInfo("Cluster pointer null!");
2237  return;
2238  }
2239 
2240  if (cluster->GetM02() != 0)
2241  fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
2242 
2243  Float_t pidlist[AliPID::kSPECIESCN+1];
2244  for (Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
2245 
2246  cluster->SetPID(pidlist);
2247 }
2248 
2269 //___________________________________________________________________________________________________________________
2271  AliVCaloCells* cells, AliVCluster * cluster,
2272  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2273  Float_t & enAfterCuts, Float_t & l0, Float_t & l1,
2274  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2275  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2276 {
2277  if (!cluster)
2278  {
2279  AliInfo("Cluster pointer null!");
2280  return;
2281  }
2282 
2283  Double_t eCell = 0.;
2284  Double_t tCell = 0.;
2285  Float_t fraction = 1.;
2286  Float_t recalFactor = 1.;
2287  Bool_t isLowGain = kFALSE;
2288 
2289  Int_t iSupMod = -1;
2290  Int_t iTower = -1;
2291  Int_t iIphi = -1;
2292  Int_t iIeta = -1;
2293  Int_t iphi = -1;
2294  Int_t ieta = -1;
2295  Double_t etai = -1.;
2296  Double_t phii = -1.;
2297 
2298  Int_t nstat = 0 ;
2299  Float_t wtot = 0.;
2300  Double_t w = 0.;
2301  Double_t etaMean = 0.;
2302  Double_t phiMean = 0.;
2303 
2304  Double_t pGlobal[3];
2305 
2306  // Loop on cells, calculate the cluster energy, in case a cut on cell energy is added,
2307  // or the non linearity correction was applied
2308  // and to check if the cluster is between 2 SM in eta
2309  Int_t iSM0 = -1;
2310  Bool_t shared = kFALSE;
2311  Float_t energy = 0;
2312 
2313  enAfterCuts = 0;
2314  l0 = 0; l1 = 0;
2315  disp = 0; dEta = 0; dPhi = 0;
2316  sEta = 0; sPhi = 0; sEtaPhi = 0;
2317 
2318  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2319  {
2320  // Get from the absid the supermodule, tower and eta/phi numbers
2321 
2322  Int_t absId = cluster->GetCellAbsId(iDigit);
2323 
2324  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2325  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2326 
2327  // Check if there are cells of different SM
2328  if (iDigit == 0 ) iSM0 = iSupMod;
2329  else if (iSupMod!= iSM0) shared = kTRUE;
2330 
2331  // Get the cell energy, if recalibration is on, apply factors
2332  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2333  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2334 
2335  if (IsRecalibrationOn())
2336  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2337 
2338  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2339  tCell = cells->GetCellTime (absId);
2340  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2341 
2342  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2343  tCell*=1e9;
2344  tCell-=fConstantTimeShift;
2345 
2346  if(eCell > 0.05) // at least a minimum 50 MeV cut
2347  energy += eCell;
2348 
2349  if(eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut)
2350  enAfterCuts += eCell;
2351  } // cell loop
2352 
2353 
2354  // Loop on cells to calculate weights and shower shape terms parameters
2355  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2356  {
2357  // Get from the absid the supermodule, tower and eta/phi numbers
2358  Int_t absId = cluster->GetCellAbsId(iDigit);
2359 
2360  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2361  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2362 
2363  //Get the cell energy, if recalibration is on, apply factors
2364  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2365  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2366 
2368  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2369 
2370  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2371  tCell = cells->GetCellTime (absId);
2372  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2373 
2374  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2375  tCell*=1e9;
2376  tCell-=fConstantTimeShift;
2377 
2378  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2379  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2380  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2381 
2382  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2383  {
2384  w = GetCellWeight(eCell, energy);
2385 
2386  // Cell index
2387  if ( fShowerShapeCellLocationType == 0 )
2388  {
2389  etai=(Double_t)ieta;
2390  phii=(Double_t)iphi;
2391  }
2392  // Cell angle location
2393  else if( fShowerShapeCellLocationType == 1 )
2394  {
2395  geom->EtaPhiFromIndex(absId, etai, phii);
2396  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2397  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2398  }
2399  else
2400  {
2401  geom->GetGlobal(absId,pGlobal);
2402 
2403  // Cell x-z location
2404  if( fShowerShapeCellLocationType == 2 )
2405  {
2406  etai = pGlobal[2];
2407  phii = pGlobal[0];
2408  }
2409  // Cell r-z location
2410  else
2411  {
2412  etai = pGlobal[2];
2413  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2414  }
2415  }
2416 
2417  if (w > 0.0)
2418  {
2419  wtot += w ;
2420  nstat++;
2421 
2422  // Shower shape
2423  sEta += w * etai * etai ;
2424  etaMean += w * etai ;
2425  sPhi += w * phii * phii ;
2426  phiMean += w * phii ;
2427  sEtaPhi += w * etai * phii ;
2428  }
2429  }
2430  else if(eCell > 0.05)
2431  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2432  } // cell loop
2433 
2434  //printf("sEta %f sPhi %f etaMean %f phiMean %f sEtaPhi %f wtot %f\n",sEta,sPhi,etaMean,phiMean,sEtaPhi, wtot);
2435 
2436  // Normalize to the weight
2437  if (wtot > 0)
2438  {
2439  etaMean /= wtot ;
2440  phiMean /= wtot ;
2441  }
2442  else
2443  AliDebug(2,Form("Wrong weight %f\n", wtot));
2444 
2445  // Loop on cells to calculate dispersion
2446  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2447  {
2448  // Get from the absid the supermodule, tower and eta/phi numbers
2449  Int_t absId = cluster->GetCellAbsId(iDigit);
2450 
2451  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2452  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2453 
2454  //Get the cell energy, if recalibration is on, apply factors
2455  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2456  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2457  if (IsRecalibrationOn())
2458  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2459 
2460  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2461  tCell = cells->GetCellTime (absId);
2462  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2463 
2464  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2465  tCell*=1e9;
2466  tCell-=fConstantTimeShift;
2467 
2468  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2469  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2470  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2471 
2472  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2473  {
2474  w = GetCellWeight(eCell,cluster->E());
2475 
2476  // Cell index
2477  if ( fShowerShapeCellLocationType == 0 )
2478  {
2479  etai=(Double_t)ieta;
2480  phii=(Double_t)iphi;
2481  }
2482  // Cell angle location
2483  else if( fShowerShapeCellLocationType == 1 )
2484  {
2485  geom->EtaPhiFromIndex(absId, etai, phii);
2486  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2487  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2488  }
2489  else
2490  {
2491  geom->GetGlobal(absId,pGlobal);
2492 
2493  // Cell x-z location
2494  if( fShowerShapeCellLocationType == 2 )
2495  {
2496  etai = pGlobal[2];
2497  phii = pGlobal[0];
2498  }
2499  // Cell r-z location
2500  else
2501  {
2502  etai = pGlobal[2];
2503  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2504  }
2505  }
2506 
2507  if (w > 0.0)
2508  {
2509  disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
2510  dEta += w * (etai-etaMean)*(etai-etaMean) ;
2511  dPhi += w * (phii-phiMean)*(phii-phiMean) ;
2512  }
2513  }
2514  else
2515  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2516  }// cell loop
2517 
2518  // Normalize to the weigth and set shower shape parameters
2519  if (wtot > 0 && nstat > 1)
2520  {
2521  disp /= wtot ;
2522  dEta /= wtot ;
2523  dPhi /= wtot ;
2524  sEta /= wtot ;
2525  sPhi /= wtot ;
2526  sEtaPhi /= wtot ;
2527 
2528  sEta -= etaMean * etaMean ;
2529  sPhi -= phiMean * phiMean ;
2530  sEtaPhi -= etaMean * phiMean ;
2531 
2532  l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2533  l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2534 
2535  //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);
2536 
2537  }
2538  else
2539  {
2540  l0 = 0. ;
2541  l1 = 0. ;
2542  dEta = 0. ; dPhi = 0. ; disp = 0. ;
2543  sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
2544  }
2545 }
2546 
2565 //___________________________________________________________________________________________________________________
2567  AliVCaloCells* cells, AliVCluster * cluster,
2568  Float_t & l0, Float_t & l1,
2569  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2570  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2571 {
2572  Float_t newEnergy = 0;
2573  Float_t cellEmin = 0.05; // 50 MeV
2574  Float_t cellTimeWindow = 1000; // open cut
2575  Int_t bc = 0;
2577  cellEmin, cellTimeWindow, bc,
2578  newEnergy, l0, l1, disp,
2579  dEta, dPhi, sEta, sPhi, sEtaPhi);
2580 }
2581 
2590 //____________________________________________________________________________________________
2592  AliVCaloCells* cells,
2593  AliVCluster * cluster)
2594 {
2595  Float_t l0 = 0., l1 = 0.;
2596  Float_t disp = 0., dEta = 0., dPhi = 0.;
2597  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2598 
2599  AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
2600  dEta, dPhi, sEta, sPhi, sEtaPhi);
2601 
2602  cluster->SetM02(l0);
2603  cluster->SetM20(l1);
2604  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2605 
2606 }
2607 
2620 //____________________________________________________________________________________________
2622  AliVCaloCells* cells, AliVCluster * cluster,
2623  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2624  Float_t & enAfterCuts)
2625 {
2626  Float_t l0 = 0., l1 = 0.;
2627  Float_t disp = 0., dEta = 0., dPhi = 0.;
2628  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2629 
2631  cellEcut, cellTimeCut, bc,
2632  enAfterCuts, l0, l1, disp,
2633  dEta, dPhi, sEta, sPhi, sEtaPhi);
2634 
2635  cluster->SetM02(l0);
2636  cluster->SetM20(l1);
2637  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2638 
2639 }
2640 
2654 //____________________________________________________________________________
2655 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
2656  TObjArray * clusterArr,
2657  const AliEMCALGeometry *geom,
2658  AliMCEvent * mc)
2659 {
2660  fMatchedTrackIndex ->Reset();
2661  fMatchedClusterIndex->Reset();
2662  fResidualPhi->Reset();
2663  fResidualEta->Reset();
2664 
2665  fMatchedTrackIndex ->Set(1000);
2666  fMatchedClusterIndex->Set(1000);
2667  fResidualPhi->Set(1000);
2668  fResidualEta->Set(1000);
2669 
2670  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
2671  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
2672 
2673  // Init the magnetic field if not already on
2674  if (!TGeoGlobalMagField::Instance()->GetField())
2675  {
2676  if (!event->InitMagneticField())
2677  {
2678  AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
2679  }
2680  } // Init mag field
2681 
2682  if (esdevent)
2683  {
2684  UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
2685  UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
2686  Bool_t desc1 = (mask1 >> 3) & 0x1;
2687  Bool_t desc2 = (mask2 >> 3) & 0x1;
2688  if (desc1==0 || desc2==0) {
2689  // AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
2690  // mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
2691  // mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
2692  fITSTrackSA=kTRUE;
2693  }
2694  }
2695 
2696  TObjArray *clusterArray = 0x0;
2697  if (!clusterArr)
2698  {
2699  clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
2700  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2701  {
2702  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2703  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells()))
2704  continue;
2705  clusterArray->AddAt(cluster,icl);
2706  }
2707  }
2708 
2709  Int_t matched=0;
2710  Double_t cv[21];
2711  TString genName;
2712  for (Int_t i=0; i<21;i++) cv[i]=0;
2713  for (Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
2714  {
2715  AliExternalTrackParam *trackParam = 0;
2716  Int_t mcLabel = -1;
2717  // If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
2718  AliESDtrack *esdTrack = 0;
2719  AliAODTrack *aodTrack = 0;
2720  if (esdevent)
2721  {
2722  esdTrack = esdevent->GetTrack(itr);
2723  if (!esdTrack) continue;
2724  if (!IsAccepted(esdTrack)) continue;
2725  if (esdTrack->Pt()<fCutMinTrackPt) continue;
2726 
2727  if ( TMath::Abs(esdTrack->Eta()) > 0.9 ) continue;
2728 
2729  // Save some time and memory in case of no DCal present
2730  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2731  {
2732  Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
2733  if ( phi <= 10 || phi >= 250 ) continue;
2734  }
2735 
2736  if (!fITSTrackSA) // if TPC Available
2737  {
2738  if ( fUseOuterTrackParam )
2739  trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetOuterParam());
2740  else
2741  trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
2742  }
2743  else
2744  trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
2745 
2746  mcLabel = TMath::Abs(esdTrack->GetLabel());
2747  }
2748 
2749  // If the input event is AOD, the starting point for extrapolation is at vertex
2750  // AOD tracks are selected according to its filterbit.
2751  else if (aodevent)
2752  {
2753  aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(itr));
2754 
2755  if (!aodTrack) AliFatal("Not a standard AOD");
2756 
2757  if (!aodTrack) continue;
2758 
2759  if (fAODTPCOnlyTracks)
2760  { // Match with TPC only tracks, default from May 2013, before filter bit 32
2761  //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
2762  if (!aodTrack->IsTPCConstrained()) continue ;
2763  }
2764  else if (fAODHybridTracks)
2765  { // Match with hybrid tracks
2766  //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
2767  if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
2768  } else
2769  { // Match with tracks on a mask
2770  //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask));
2771  if (!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks
2772  }
2773 
2774  if (aodTrack->Pt() < fCutMinTrackPt) continue;
2775 
2776  if ( TMath::Abs(aodTrack->Eta()) > 0.9 ) continue;
2777 
2778  // Save some time and memory in case of no DCal present
2779  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2780  {
2781  Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
2782  if ( phi <= 10 || phi >= 250 ) continue;
2783  }
2784 
2785  Double_t pos[3],mom[3];
2786  aodTrack->GetXYZ(pos);
2787  aodTrack->GetPxPyPz(mom);
2788  AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",
2789  itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
2790 
2791  trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
2792 
2793  mcLabel = TMath::Abs(aodTrack->GetLabel());
2794  }
2795 
2796  //Return if the input data is not "AOD" or "ESD"
2797  else
2798  {
2799  AliWarning("Wrong input data type! Should be \"AOD\" or \"ESD\" ");
2800  if (clusterArray)
2801  {
2802  clusterArray->Clear();
2803  delete clusterArray;
2804  }
2805  return;
2806  }
2807 
2808  if (!trackParam) continue;
2809 
2810  //
2811  // Check if track comes from a particular MC generator, do not include it
2812  // if it is not a selected one
2813  //
2814  if( mc && fMCGenerToAcceptForTrack && fNMCGenerToAccept > 0 )
2815  {
2816  mc->GetCocktailGenerator(mcLabel,genName);
2817 
2818  Bool_t generOK = kFALSE;
2819  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
2820  {
2821  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
2822  }
2823 
2824  if ( !generOK ) continue;
2825  }
2826 
2827  // Extrapolate the track to EMCal surface, see AliEMCALRecoUtilsBase
2828  AliExternalTrackParam emcalParam(*trackParam);
2829  Float_t eta, phi, pt;
2830  if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2831  {
2832  if (aodevent && trackParam) delete trackParam;
2833  if (fITSTrackSA && trackParam) delete trackParam;
2834  continue;
2835  }
2836 
2837  if ( TMath::Abs(eta) > 0.75 )
2838  {
2839  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2840  continue;
2841  }
2842 
2843  // Save some time and memory in case of no DCal present
2844  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2845  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2846  {
2847  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2848  continue;
2849  }
2850 
2851  //Find matched clusters
2852  Int_t index = -1;
2853  Float_t dEta = -999, dPhi = -999;
2854  if (!clusterArr)
2855  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
2856  else
2857  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr , dEta, dPhi);
2858 
2859 
2860  if (index>-1)
2861  {
2862  fMatchedTrackIndex ->AddAt(itr,matched);
2863  fMatchedClusterIndex ->AddAt(index,matched);
2864  fResidualEta ->AddAt(dEta,matched);
2865  fResidualPhi ->AddAt(dPhi,matched);
2866  matched++;
2867  }
2868 
2869  if (aodevent && trackParam) delete trackParam;
2870  if (fITSTrackSA && trackParam) delete trackParam;
2871  }//track loop
2872 
2873  if (clusterArray)
2874  {
2875  clusterArray->Clear();
2876  delete clusterArray;
2877  }
2878 
2879  AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
2880 
2881  fMatchedTrackIndex ->Set(matched);
2882  fMatchedClusterIndex ->Set(matched);
2883  fResidualPhi ->Set(matched);
2884  fResidualEta ->Set(matched);
2885 }
2886 
2898 //________________________________________________________________________________
2900  const AliVEvent *event,
2901  const AliEMCALGeometry *geom,
2902  Float_t &dEta, Float_t &dPhi)
2903 {
2904  Int_t index = -1;
2905 
2906  if ( TMath::Abs(track->Eta()) > 0.9 ) return index;
2907 
2908  // Save some time and memory in case of no DCal present
2909  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2910  {
2911  Double_t phiV = track->Phi()*TMath::RadToDeg();
2912  if ( phiV <= 10 || phiV >= 250 ) return index;
2913  }
2914 
2915  AliExternalTrackParam *trackParam = 0;
2916  if (!fITSTrackSA) // If TPC
2917  {
2918  if ( fUseOuterTrackParam )
2919  trackParam = const_cast<AliExternalTrackParam*>(track->GetOuterParam());
2920  else
2921  trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
2922  }
2923  else
2924  trackParam = new AliExternalTrackParam(*track);
2925 
2926  if (!trackParam) return index;
2927  AliExternalTrackParam emcalParam(*trackParam);
2928 
2929  Float_t eta, phi, pt;
2930  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2931  {
2932  if (fITSTrackSA) delete trackParam;
2933  return index;
2934  }
2935 
2936  if ( TMath::Abs(eta) > 0.75 )
2937  {
2938  if (fITSTrackSA) delete trackParam;
2939  return index;
2940  }
2941 
2942  // Save some time and memory in case of no DCal present
2943  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2944  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2945  {
2946  if (fITSTrackSA) delete trackParam;
2947  return index;
2948  }
2949 
2950  TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
2951 
2952  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2953  {
2954  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2955  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
2956  clusterArr->AddAt(cluster,icl);
2957  }
2958 
2959  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
2960  clusterArr->Clear();
2961  delete clusterArr;
2962  if (fITSTrackSA) delete trackParam;
2963 
2964  return index;
2965 }
2966 
2977 //_______________________________________________________________________________________________
2978 Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
2979  AliExternalTrackParam *trkParam,
2980  const TObjArray * clusterArr,
2981  Float_t &dEta, Float_t &dPhi)
2982 {
2983  dEta=-999, dPhi=-999;
2984  Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
2985  Int_t index = -1;
2986  Float_t tmpEta=-999, tmpPhi=-999;
2987 
2988  Double_t exPos[3] = {0.,0.,0.};
2989  if (!emcalParam->GetXYZ(exPos)) return index;
2990 
2991  Float_t clsPos[3] = {0.,0.,0.};
2992  for (Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
2993  {
2994  AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
2995 
2996  if (!cluster || !cluster->IsEMCAL()) continue;
2997 
2998  cluster->GetPosition(clsPos);
2999 
3000  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));
3001  if (dR > fClusterWindow) continue;
3002 
3003  AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
3004 
3005  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
3006 
3007  if (fCutEtaPhiSum)
3008  {
3009  Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
3010  if (tmpR<dRMax)
3011  {
3012  dRMax=tmpR;
3013  dEtaMax=tmpEta;
3014  dPhiMax=tmpPhi;
3015  index=icl;
3016  }
3017  }
3018  else if (fCutEtaPhiSeparate)
3019  {
3020  if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
3021  {
3022  dEtaMax = tmpEta;
3023  dPhiMax = tmpPhi;
3024  index=icl;
3025  }
3026  }
3027  else
3028  {
3029  AliError("Please specify your cut criteria");
3030  AliError("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()");
3031  AliError("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()");
3032  return index;
3033  }
3034  }
3035 
3036  dEta=dEtaMax;
3037  dPhi=dPhiMax;
3038 
3039  return index;
3040 }
3041 
3054 //---------------------------------------------------------------------------------
3055 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
3056  const AliVCluster *cluster,
3057  Float_t &tmpEta,
3058  Float_t &tmpPhi)
3059 {
3060  return AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
3061 }
3062 
3072 //_______________________________________________________________________
3074  Float_t &dEta, Float_t &dPhi)
3075 {
3076  if (FindMatchedPosForCluster(clsIndex) >= 999)
3077  {
3078  AliDebug(2,"No matched tracks found!\n");
3079  dEta=999.;
3080  dPhi=999.;
3081  return;
3082  }
3083 
3084  dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
3085  dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
3086 }
3087 
3097 //______________________________________________________________________________________________
3099 {
3100  if (FindMatchedPosForTrack(trkIndex) >= 999)
3101  {
3102  AliDebug(2,"No matched cluster found!\n");
3103  dEta=999.;
3104  dPhi=999.;
3105  return;
3106  }
3107 
3108  dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
3109  dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
3110 }
3111 
3118 //__________________________________________________________
3120 {
3121  if (IsClusterMatched(clsIndex))
3122  return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
3123  else
3124  return -1;
3125 }
3126 
3133 //__________________________________________________________
3135 {
3136  if (IsTrackMatched(trkIndex))
3137  return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
3138  else
3139  return -1;
3140 }
3141 
3149 //______________________________________________________________
3151 {
3152  if (FindMatchedPosForCluster(clsIndex) < 999)
3153  return kTRUE;
3154  else
3155  return kFALSE;
3156 }
3157 
3165 //____________________________________________________________
3167 {
3168  if (FindMatchedPosForTrack(trkIndex) < 999)
3169  return kTRUE;
3170  else
3171  return kFALSE;
3172 }
3173 
3181 //______________________________________________________________________
3183 {
3184  Float_t tmpR = fCutR;
3185  UInt_t pos = 999;
3186 
3187  for (Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
3188  {
3189  if (fMatchedClusterIndex->At(i)==clsIndex)
3190  {
3191  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3192  if (r<tmpR)
3193  {
3194  pos=i;
3195  tmpR=r;
3196  AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3197  fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3198  }
3199  }
3200  }
3201  return pos;
3202 }
3203 
3211 //____________________________________________________________________
3213 {
3214  Float_t tmpR = fCutR;
3215  UInt_t pos = 999;
3216 
3217  for (Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
3218  {
3219  if (fMatchedTrackIndex->At(i)==trkIndex)
3220  {
3221  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3222  if (r<tmpR)
3223  {
3224  pos=i;
3225  tmpR=r;
3226  AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3227  fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3228  }
3229  }
3230  }
3231  return pos;
3232 }
3233 
3248 //__________________________________________________________________________
3250  const AliEMCALGeometry *geom,
3251  AliVCaloCells* cells, Int_t bc)
3252 {
3253  Bool_t isGood=kTRUE;
3254 
3255  if (!cluster || !cluster->IsEMCAL()) return kFALSE;
3256  if (ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
3257  if (!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
3258  if (IsExoticCluster(cluster, cells,bc)) return kFALSE;
3259 
3260  return isGood;
3261 }
3262 
3273 //__________________________________________________________
3275 {
3276  UInt_t status = esdTrack->GetStatus();
3277 
3278  Int_t nClustersITS = esdTrack->GetITSclusters(0);
3279  Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
3280 
3281  Float_t chi2PerClusterITS = -1;
3282  Float_t chi2PerClusterTPC = -1;
3283  if (nClustersITS!=0)
3284  chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
3285  if (nClustersTPC!=0)
3286  chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
3287 
3288  //
3289  // DCA cuts
3290  // Only to be used for primary
3291  //
3292  if ( fTrackCutsType == kGlobalCut )
3293  {
3294  Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01);
3295  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
3296 
3297  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3298 
3299  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3300  }
3301  else if( fTrackCutsType == kGlobalCut2011 )
3302  {
3303  Float_t maxDCAToVertexXYPtDep = 0.0105 + 0.0350/TMath::Power(esdTrack->Pt(),1.1);
3304  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2011()
3305 
3306  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3307 
3308  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3309  }
3310 
3311  Float_t b[2];
3312  Float_t bCov[3];
3313  esdTrack->GetImpactParameters(b,bCov);
3314  if (bCov[0]<=0 || bCov[2]<=0)
3315  {
3316  AliDebug(1, "Estimated b resolution lower or equal zero!");
3317  bCov[0]=0; bCov[2]=0;
3318  }
3319 
3320  Float_t dcaToVertexXY = b[0];
3321  Float_t dcaToVertexZ = b[1];
3322  Float_t dcaToVertex = -1;
3323 
3324  if (fCutDCAToVertex2D)
3325  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY / fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY +
3326  dcaToVertexZ*dcaToVertexZ / fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ );
3327  else
3328  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
3329 
3330  // cut the track?
3331 
3332  Bool_t cuts[kNCuts];
3333  for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
3334 
3335  // track quality cuts
3336  if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
3337  cuts[0]=kTRUE;
3338  if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
3339  cuts[1]=kTRUE;
3340  if (nClustersTPC<fCutMinNClusterTPC)
3341  cuts[2]=kTRUE;
3342  if (nClustersITS<fCutMinNClusterITS)
3343  cuts[3]=kTRUE;
3344  if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
3345  cuts[4]=kTRUE;
3346  if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
3347  cuts[5]=kTRUE;
3348  if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
3349  cuts[6]=kTRUE;
3350  if (fCutDCAToVertex2D && dcaToVertex > 1)
3351  cuts[7] = kTRUE;
3352  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
3353  cuts[8] = kTRUE;
3354  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
3355  cuts[9] = kTRUE;
3356 
3358  {
3359  //Require at least one SPD point + anything else in ITS
3360  if ( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
3361  cuts[10] = kTRUE;
3362  }
3363 
3364  // ITS
3366  {
3367  if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin))
3368  {
3369  // TPC tracks
3370  cuts[11] = kTRUE;
3371  }
3372  else
3373  {
3374  // ITS standalone tracks
3376  {
3377  if (status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
3378  }
3379  else if (fCutRequireITSpureSA)
3380  {
3381  if (!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
3382  }
3383  }
3384  }
3385 
3386  Bool_t cut=kFALSE;
3387  for (Int_t i=0; i<kNCuts; i++)
3388  if (cuts[i]) { cut = kTRUE ; }
3389 
3390  // cut the track
3391  if (cut)
3392  return kFALSE;
3393  else
3394  return kTRUE;
3395 }
3396 
3402 //_____________________________________
3404 {
3405  switch (fTrackCutsType)
3406  {
3407  case kTPCOnlyCut:
3408  {
3409  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()"));
3410  //TPC
3411  SetMinNClustersTPC(70);
3413  SetAcceptKinkDaughters(kFALSE);
3414  SetRequireTPCRefit(kFALSE);
3415 
3416  //ITS
3417  SetRequireITSRefit(kFALSE);
3418  SetMaxDCAToVertexZ(3.2);
3419  SetMaxDCAToVertexXY(2.4);
3420  SetDCAToVertex2D(kTRUE);
3421 
3422  break;
3423  }
3424 
3425  case kGlobalCut:
3426  {
3427  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE)"));
3428  //TPC
3429  SetMinNClustersTPC(70);
3431  SetAcceptKinkDaughters(kFALSE);
3432  SetRequireTPCRefit(kTRUE);
3433 
3434  //ITS
3435  SetRequireITSRefit(kTRUE);
3436  SetMaxDCAToVertexZ(2);
3438  SetDCAToVertex2D(kFALSE);
3439 
3440  break;
3441  }
3442 
3443  case kLooseCut:
3444  {
3445  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
3446  SetMinNClustersTPC(50);
3447  SetAcceptKinkDaughters(kTRUE);
3448 
3449  break;
3450  }
3451 
3452  case kITSStandAlone:
3453  {
3454  AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
3455  SetRequireITSRefit(kTRUE);
3456  SetRequireITSStandAlone(kTRUE);
3457  SetITSTrackSA(kTRUE);
3458  break;
3459  }
3460 
3461  case kGlobalCut2011:
3462  {
3463  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE)"));
3464  //TPC
3465  SetMinNClustersTPC(50);
3467  SetAcceptKinkDaughters(kFALSE);
3468  SetRequireTPCRefit(kTRUE);
3469 
3470  //ITS
3471  SetRequireITSRefit(kTRUE);
3472  SetMaxDCAToVertexZ(2);
3474  SetDCAToVertex2D(kFALSE);
3475 
3476  break;
3477  }
3478 
3479  case kLooseCutWithITSrefit:
3480  {
3481  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut plus ITSrefit"));
3482  SetMinNClustersTPC(50);
3483  SetAcceptKinkDaughters(kTRUE);
3484  SetRequireITSRefit(kTRUE);
3485 
3486  break;
3487  }
3488  }
3489 }
3490 
3496 //________________________________________________________________________
3498 {
3499  Int_t nTracks = event->GetNumberOfTracks();
3500  for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
3501  {
3502  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
3503  if (!track)
3504  {
3505  AliWarning(Form("Could not receive track %d", iTrack));
3506  continue;
3507  }
3508 
3509  Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
3510  track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
3511 
3512  /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
3513  AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
3514  if (esdtrack)
3515  {
3516  if (matchClusIndex != -1)
3517  esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
3518  else
3519  esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3520  }
3521  else
3522  {
3523  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
3524  if (matchClusIndex != -1)
3525  aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
3526  else
3527  aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3528  }
3529  }
3530  AliDebug(2,"Track matched to closest cluster");
3531 }
3532 
3539 //_________________________________________________________________________
3541 {
3542  for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
3543  {
3544  AliVCluster *cluster = event->GetCaloCluster(iClus);
3545  if (!cluster->IsEMCAL())
3546  continue;
3547 
3548  //
3549  // Remove old matches in cluster
3550  //
3551  if(cluster->GetNTracksMatched() > 0)
3552  {
3553  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
3554  {
3555  TArrayI arrayTrackMatched(0);
3556  ((AliESDCaloCluster*)cluster)->AddTracksMatched(arrayTrackMatched);
3557  }
3558  else
3559  {
3560  for(Int_t iTrack = 0; iTrack < cluster->GetNTracksMatched(); iTrack++)
3561  {
3562  ((AliAODCaloCluster*)cluster)->RemoveTrackMatched((TObject*)((AliAODCaloCluster*)cluster)->GetTrackMatched(iTrack));
3563  }
3564  }
3565  }
3566 
3567  //
3568  // Find new matches and put them in the cluster
3569  //
3570  Int_t nTracks = event->GetNumberOfTracks();
3571  TArrayI arrayTrackMatched(nTracks);
3572 
3573  // Get the closest track matched to the cluster
3574  Int_t nMatched = 0;
3575  Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
3576  if (matchTrackIndex != -1)
3577  {
3578  arrayTrackMatched[nMatched] = matchTrackIndex;
3579  nMatched++;
3580  }
3581 
3582  // Get all other tracks matched to the cluster
3583  for (Int_t iTrk=0; iTrk<nTracks; ++iTrk)
3584  {
3585  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
3586 
3587  if( !track ) continue;
3588 
3589  if ( iTrk == matchTrackIndex ) continue;
3590 
3591  if ( track->GetEMCALcluster() == iClus )
3592  {
3593  arrayTrackMatched[nMatched] = iTrk;
3594  ++nMatched;
3595  }
3596  }
3597 
3598  AliDebug(2,Form("cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]));
3599 
3600  arrayTrackMatched.Set(nMatched);
3601  AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
3602  if (esdcluster)
3603  esdcluster->AddTracksMatched(arrayTrackMatched);
3604  else if ( nMatched > 0 )
3605  {
3606  AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
3607  if (aodcluster)
3608  {
3609  aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
3610  //AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(event->GetTrack(arrayTrackMatched.At(0)));
3611  //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());
3612  //printf("With specs: pt %.4f, eta %.4f, phi %.4f\n",aodtrack->Pt(),aodtrack->Eta(), aodtrack->Phi());
3613  }
3614  }
3615 
3616  Float_t eta= -999, phi = -999;
3617  if (matchTrackIndex != -1)
3618  GetMatchedResiduals(iClus, eta, phi);
3619 
3620  cluster->SetTrackDistance(phi, eta);
3621  }
3622 
3623  AliDebug(2,"Cluster matched to tracks");
3624 }
3625 
3628  else {
3629  fEMCALRecalibrationFactors = new TObjArray(map->GetEntries());
3630  fEMCALRecalibrationFactors->SetOwner(true);
3631  }
3632  if(!fEMCALRecalibrationFactors->IsOwner()){
3633  // Must claim ownership since the new objects are owend by this instance
3634  fEMCALRecalibrationFactors->SetOwner(kTRUE);
3635  }
3636  for(int i = 0; i < map->GetEntries(); i++){
3637  TH2F *hist = dynamic_cast<TH2F *>(map->At(i));
3638  if(!hist) continue;
3639  this->SetEMCALChannelRecalibrationFactors(i, hist);
3640  }
3641 }
3642 
3646  fEMCALRecalibrationFactors->SetOwner(true);
3647  }
3648  if(fEMCALRecalibrationFactors->GetEntries() <= iSM) fEMCALRecalibrationFactors->Expand(iSM+1);
3649  if(fEMCALRecalibrationFactors->At(iSM)) fEMCALRecalibrationFactors->RemoveAt(iSM);
3650  TH2F *clone = new TH2F(*h);
3651  clone->SetDirectory(NULL);
3652  fEMCALRecalibrationFactors->AddAt(clone,iSM);
3653 }
3654 
3657  else {
3658  fEMCALBadChannelMap = new TObjArray(map->GetEntries());
3659  fEMCALBadChannelMap->SetOwner(true);
3660  }
3661  if(!fEMCALBadChannelMap->IsOwner()) {
3662  // Must claim ownership since the new objects are owend by this instance
3663  fEMCALBadChannelMap->SetOwner(true);
3664  }
3665  for(int i = 0; i < map->GetEntries(); i++){
3666  TH2I *hist = dynamic_cast<TH2I *>(map->At(i));
3667  if(!hist) continue;
3668  this->SetEMCALChannelStatusMap(i, hist);
3669  }
3670 }
3671 
3673  if(!fEMCALBadChannelMap){
3674  fEMCALBadChannelMap = new TObjArray(iSM);
3675  fEMCALBadChannelMap->SetOwner(true);
3676  }
3677  if(fEMCALBadChannelMap->GetEntries() <= iSM) fEMCALBadChannelMap->Expand(iSM+1);
3678  if(fEMCALBadChannelMap->At(iSM)) fEMCALBadChannelMap->RemoveAt(iSM);
3679  TH2I *clone = new TH2I(*h);
3680  clone->SetDirectory(NULL);
3681  fEMCALBadChannelMap->AddAt(clone,iSM);
3682 }
3683 
3686  else {
3687  fEMCALTimeRecalibrationFactors = new TObjArray(map->GetEntries());
3688  fEMCALTimeRecalibrationFactors->SetOwner(true);
3689  }
3690  if(!fEMCALTimeRecalibrationFactors->IsOwner()) {
3691  // Must claim ownership since the new objects are owend by this instance
3692  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
3693  }
3694  for(int i = 0; i < map->GetEntries(); i++){
3695  TH1F *hist = dynamic_cast<TH1F *>(map->At(i));
3696  if(!hist) continue;
3698  }
3699 }
3700 
3704  fEMCALTimeRecalibrationFactors->SetOwner(true);
3705  }
3706  if(fEMCALTimeRecalibrationFactors->GetEntries() <= bc) fEMCALTimeRecalibrationFactors->Expand(bc+1);
3708  TH1F *clone = new TH1F(*h);
3709  clone->SetDirectory(NULL);
3710  fEMCALTimeRecalibrationFactors->AddAt(clone,bc);
3711 }
3712 
3715  else {
3716  fEMCALL1PhaseInTimeRecalibration = new TObjArray(map->GetEntries());
3717  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3718  }
3719  if(!fEMCALL1PhaseInTimeRecalibration->IsOwner()) {
3720  // Must claim ownership since the new objects are owend by this instance
3721  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3722  }
3723  for(int i = 0; i < map->GetEntries(); i++) {
3724  TH1C *hist = dynamic_cast<TH1C *>(map->At(i));
3725  if(!hist) continue;
3727  }
3728 }
3729 
3733  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3734  }
3736  TH1C *clone = new TH1C(*h);
3737  clone->SetDirectory(NULL);
3738  fEMCALL1PhaseInTimeRecalibration->AddAt(clone,0);
3739 }
3740 
3744 //___________________________________________________
3746 {
3747  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3748  printf("AliEMCALRecoUtils Settings: \n");
3749  printf("\tMisalignment shifts\n");
3750  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,
3752  fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
3753  printf("\tNon linearity function %d, parameters:\n", fNonLinearityFunction);
3754  if (fNonLinearityFunction != 3) // print only if not kNoCorrection
3755  for (Int_t i=0; i<10; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
3756 
3757  printf("\tPosition Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
3758 
3759  printf("\tMatching criteria: ");
3760  if (fCutEtaPhiSum) {
3761  printf("\t\tsqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
3762  } else if (fCutEtaPhiSeparate) {
3763  printf("\t\tdEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
3764  } else {
3765  printf("\t\tError\n");
3766  printf("\t\tplease specify your cut criteria\n");
3767  printf("\t\tTo cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
3768  printf("\t\tTo cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
3769  }
3770 
3771  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);
3772  printf("\tCluster selection window: dR < %2.0f\n",fClusterWindow);
3773 
3774  printf("\tTrack cuts: \n");
3775  printf("\t\tMinimum track pT: %1.2f\n",fCutMinTrackPt);
3776  printf("\t\tAOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask);
3777  printf("\t\tTPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
3778  printf("\t\tAcceptKinks = %d\n",fCutAcceptKinkDaughters);
3779  printf("\t\tMinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
3780  printf("\t\tMaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
3781  printf("\t\tDCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
3782  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3783 }
void RecalculateCellLabelsRemoveAddedGenerator(Int_t absID, AliVCluster *clus, AliMCEvent *mc, Float_t &amp, TArrayI &labeArr, TArrayF &eDepArr) const
Bool_t IsL1PhaseInTimeRecalibrationOn() const
Double_t fMass
Mass hypothesis of the track.
Bool_t fMCGenerToAcceptForTrack
Activate the removal of tracks entering the track matching that come from a particular generator...
double Double_t
Definition: External.C:58
Bool_t CheckCellFiducialRegion(const AliEMCALGeometry *geom, const AliVCluster *cluster, AliVCaloCells *cells)
Bool_t fAODTPCOnlyTracks
Match with TPC only tracks.
Bool_t fCutRequireTPCRefit
Require TPC refit.
Definition: External.C:236
void SetRequireITSStandAlone(Bool_t b=kFALSE)
Bool_t IsClusterMatched(Int_t clsIndex) const
Float_t GetEMCALChannelTimeRecalibrationFactor(Int_t bc, Int_t absID, Bool_t isLGon=kFALSE) const
TArrayI * fMatchedClusterIndex
Array that stores indexes of matched clusters.
TObjArray * fEMCALTimeRecalibrationFactors
Array of histograms with map of time recalibration factors, EMCAL.
void SetMaxDCAToVertexXY(Float_t dist=1e10)
TRandom3 fRandom
Random generator for cluster energy smearing.
Bool_t fTimeRecalibration
Switch on or off the time recalibration.
Bool_t fCutEtaPhiSum
Place cut on sqrt(dEta^2+dPhi^2)
Int_t FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam, AliExternalTrackParam *trkParam, const TObjArray *clusterArr, Float_t &dEta, Float_t &dPhi)
void InitEMCALL1PhaseInTimeRecalibration()
Int_t fCutMinNClusterTPC
Min number of tpc clusters.
Float_t fCutMaxDCAToVertexZ
Track-to-vertex cut in max absolute distance in z-plane.
void SetRequireITSRefit(Bool_t b=kFALSE)
energy
Definition: HFPtSpectrum.C:44
void SetEMCALL1PhaseInTimeRecalibrationForSM(Int_t iSM, Int_t c=0)
void SetDCAToVertex2D(Bool_t b=kFALSE)
void SetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow, Double_t c=1)
Int_t fNonLinearThreshold
Non linearity threshold value for kBeamTest non linearity function.
Bool_t fITSTrackSA
If track matching is to be done with ITS tracks standing alone, ESDs.
void SetITSTrackSA(Bool_t isITS)
Bool_t IsExoticCluster(const AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=0)
Float_t fMisalRotShift[15]
Cluster position rotation shift parameters.
void SetMinNClustersTPC(Int_t min=-1)
Bool_t fRecalDistToBadChannels
Calculate distance from highest energy tower of cluster to closes bad channel.
Bool_t fRecalibration
Switch on or off the recalibration.
TArrayI * fMatchedTrackIndex
Array that stores indexes of matched tracks.
void SetMaxDCAToVertexZ(Float_t dist=1e10)
Bool_t IsBadChannelsRemovalSwitchedOn() const
Bool_t IsExoticCell(Int_t absId, AliVCaloCells *cells, Int_t bc=-1)
void RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
Float_t fExoticCellDiffTime
If time of candidate to exotic and close cell is too different (in ns), it must be noisy...
TArrayF * fResidualPhi
Array that stores the residual phi.
void RecalibrateCells(AliVCaloCells *cells, Int_t bc)
Double_t fCutMinTrackPt
Cut on track pT.
Bool_t IsRecalibrationOn() const
Some utilities for cluster and cell treatment.
Float_t fCutMaxChi2PerClusterITS
Max its fit chi2 per its cluster.
void RecalculateClusterShowerShapeParametersWithCellCuts(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster, Float_t cellEcut, Float_t cellTimeCut, Int_t bc, Float_t &enAfterCuts)
Bool_t fRemoveBadChannels
Check the channel status provided and remove clusters with bad channels.
Bool_t IsAccepted(AliESDtrack *track)
Float_t fMisalTransShift[15]
Cluster position translation shift parameters.
void SetAcceptKinkDaughters(Bool_t b=kTRUE)
void RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
Double_t fClusterWindow
Select clusters in the window to be matched.
Bool_t fCutRequireITSRefit
Require ITS refit.
Float_t fNonLinearityParams[10]
Parameters for the non linearity function.
Float_t GetCellWeight(Float_t eCell, Float_t eCluster) const
Bool_t fUseL1PhaseInTimeRecalibration
Switch on or off the L1 phase in time recalibration.
Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow) const
void SetClusterMatchedToTrack(const AliVEvent *event)
Bool_t fCutDCAToVertex2D
If true a 2D DCA cut is made.
Int_t GetMatchedTrackIndex(Int_t clsIndex)
void RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
int Int_t
Definition: External.C:63
void Print(const Option_t *) const
UInt_t FindMatchedPosForCluster(Int_t clsIndex) const
Bool_t AcceptCalibrateCell(Int_t absId, Int_t bc, Float_t &amp, Double_t &time, AliVCaloCells *cells)
unsigned int UInt_t
Definition: External.C:33
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 SetMaxChi2PerClusterTPC(Float_t max=1e10)
Float_t fConstantTimeShift
Apply a 600 ns (+15.8) time shift in case of simulation, shift in ns.
Int_t GetMatchedClusterIndex(Int_t trkIndex)
Int_t fNonLinearityFunction
Non linearity function choice, see enum NonlinearityFunctions.
Bool_t fCutRequireITSStandAlone
Require ITSStandAlone.
TObjArray * fEMCALL1PhaseInTimeRecalibration
Histogram with map of L1 phase per SM, EMCAL.
TObjArray * fEMCALRecalibrationFactors
Array of histograms with map of recalibration factors, EMCAL.
Bool_t fUseRunCorrectionFactors
Use Run Dependent Correction.
void SetTracksMatchedToCluster(const AliVEvent *event)
Float_t fCutMaxChi2PerClusterTPC
Max tpc fit chi2 per tpc cluster.
Float_t fExoticCellMinAmplitude
Check for exotic only if amplitud is larger than this value.
void RecalibrateClusterEnergy(const AliEMCALGeometry *geom, AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=-1)
Float_t fW0
Energy weight used in cluster position and shower shape calculations.
Bool_t fAODHybridTracks
Match with hybrid.
void GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
AliEMCALRecoUtils & operator=(const AliEMCALRecoUtils &)
void SetRequireTPCRefit(Bool_t b=kFALSE)
TArrayF * fResidualEta
Array that stores the residual eta.
short Short_t
Definition: External.C:23
TString fMCGenerToAccept[5]
List with name of generators that should not be included.
Bool_t fSmearClusterEnergy
Smear cluster energy, to be done only for simulated data to match real data.
Bool_t IsGoodCluster(AliVCluster *cluster, const AliEMCALGeometry *geom, AliVCaloCells *cells, Int_t bc=-1)
void GetMaxEnergyCell(const AliEMCALGeometry *geom, AliVCaloCells *cells, const AliVCluster *clu, Int_t &absId, Int_t &iSupMod, Int_t &ieta, Int_t &iphi, Bool_t &shared)
void InitEMCALTimeRecalibrationFactors()
Double_t fStepSurface
Length of step to extrapolate tracks to EMCal surface.
UInt_t fAODFilterMask
Filter mask to select AOD tracks. Refer to $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C.
static Bool_t ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, const AliVCluster *cluster, Double_t mass, Double_t step, Float_t &tmpEta, Float_t &tmpPhi)
void GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
Float_t fSmearClusterParam[3]
Energy smearing parameters.
Bool_t fCutAcceptKinkDaughters
Accepting kink daughters?
void RecalculateClusterPosition(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
Float_t fExoticCellFraction
Good cell if fraction < 1-ecross/ecell.
void RecalculateClusterPID(AliVCluster *cluster)
Int_t fPosAlgo
Position recalculation algorithm, see enum PositionAlgorithms.
void SetEMCALL1PhaseInTimeRecalibrationForAllSM(const TObjArray *map)
AliEMCALPIDUtils * fPIDUtils
Recalculate PID parameters.
Bool_t fRejectExoticCells
Remove exotic cells.
Int_t FindMatchedClusterInEvent(const AliESDtrack *track, const AliVEvent *event, const AliEMCALGeometry *geom, Float_t &dEta, Float_t &dPhi)
Bool_t IsTrackMatched(Int_t trkIndex) const
void RecalibrateCellTimeL1Phase(Int_t iSM, Int_t bc, Double_t &time) const
Float_t fCutMaxDCAToVertexXY
Track-to-vertex cut in max absolute distance in xy-plane.
void FindMatches(AliVEvent *event, TObjArray *clusterArr=0x0, const AliEMCALGeometry *geom=0x0, AliMCEvent *mc=0x0)
unsigned short UShort_t
Definition: External.C:28
TObjArray * fEMCALBadChannelMap
Array of histograms with map of bad channels, EMCAL.
void SetEMCALChannelTimeRecalibrationFactors(const TObjArray *map)
Double_t fStepCluster
Length of step to extrapolate tracks to clusters.
Bool_t fNoEMCALBorderAtEta0
Do fiducial cut in EMCAL region eta = 0?
const char Option_t
Definition: External.C:48
Int_t fNMCGenerToAccept
Number of MC generators that should not be included in analysis.
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.
Int_t GetEMCALL1PhaseInTimeRecalibrationForSM(Int_t iSM) const
bool Bool_t
Definition: External.C:53
Int_t fCutMinNClusterITS
Min number of its clusters.
Float_t CorrectClusterEnergyLinearity(AliVCluster *clu)
UInt_t FindMatchedPosForTrack(Int_t trkIndex) const
Bool_t ClusterContainsBadChannel(const AliEMCALGeometry *geom, const UShort_t *cellList, Int_t nCells)
Double_t fEMCalSurfaceDistance
EMCal surface distance (= 430 by default, the last 10 cm are propagated on a cluster-track pair basis...
void SetEMCALChannelRecalibrationFactors(const TObjArray *map)
Float_t GetECross(Int_t absID, Double_t tcell, AliVCaloCells *cells, Int_t bc)
Bool_t fLowGain
Switch on or off calibration with low gain channels.
void SetEMCALChannelStatusMap(const TObjArray *map)
Float_t fCutPhi
dPhi cut on matching
Int_t fTrackCutsType
ESD track cuts type for matching, see enum TrackCutsType.
void RecalibrateCellTime(Int_t absId, Int_t bc, Double_t &time, Bool_t isLGon=kFALSE) const
Bool_t fCutEtaPhiSeparate
Cut on dEta and dPhi separately.
Float_t fCutEta
dEta cut on matching
Float_t SmearClusterEnergy(const AliVCluster *clu)