AliPhysics  aaf9c62 (aaf9c62)
AliEMCALRecoUtils.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // ROOT includes
17 #include <TH2F.h>
18 #include <TArrayI.h>
19 #include <TArrayF.h>
20 #include <TObjArray.h>
21 
22 // STEER includes
23 #include "AliVCluster.h"
24 #include "AliVCaloCells.h"
25 #include "AliLog.h"
26 #include "AliPID.h"
27 #include "AliESDEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliESDtrack.h"
30 #include "AliAODTrack.h"
31 #include "AliExternalTrackParam.h"
32 #include "AliESDfriendTrack.h"
33 #include "AliTrackerBase.h"
34 #include "AliMCEvent.h"
35 
36 // EMCAL includes
37 #include "AliEMCALRecoUtils.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliTrackerBase.h"
40 #include "AliEMCALPIDUtils.h"
41 
43 ClassImp(AliEMCALRecoUtils) ;
45 
51 //_____________________________________
53  fParticleType(0), fPosAlgo(0),
54  fW0(0), fShowerShapeCellLocationType(0),
55  fNonLinearityFunction(0), fNonLinearThreshold(0),
56  fSmearClusterEnergy(kFALSE), fRandom(),
57  fCellsRecalibrated(kFALSE), fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
58  fConstantTimeShift(0), fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(), fLowGain(kFALSE),
59  fUseL1PhaseInTimeRecalibration(kFALSE), fEMCALL1PhaseInTimeRecalibration(),
60  fUseRunCorrectionFactors(kFALSE),
61  fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
62  fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
63  fRejectExoticCluster(kFALSE), fRejectExoticCells(kFALSE),
64  fExoticCellFraction(0), fExoticCellDiffTime(0), fExoticCellMinAmplitude(0),
65  fPIDUtils(), fAODFilterMask(0),
66  fAODHybridTracks(0), fAODTPCOnlyTracks(0),
67  fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
68  fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kFALSE), fCutEtaPhiSeparate(kFALSE),
69  fCutR(0), fCutEta(0), fCutPhi(0),
70  fClusterWindow(0), fMass(0),
71  fStepSurface(0), fStepCluster(0),
72  fITSTrackSA(kFALSE), fEMCalSurfaceDistance(440.),
73  fTrackCutsType(0), fCutMinTrackPt(0), fCutMinNClusterTPC(0),
74  fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
75  fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
76  fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0), fCutDCAToVertex2D(kFALSE),
77  fCutRequireITSStandAlone(kFALSE), fCutRequireITSpureSA(kFALSE),
78  fNMCGenerToAccept(0), fMCGenerToAcceptForTrack(1)
79 {
80  // Init parameters
82 
83  for(Int_t j = 0; j < 5; j++) fMCGenerToAccept[j] = "";
84 
85  // Track matching arrays init
88  fResidualPhi = new TArrayF();
89  fResidualEta = new TArrayF();
90  fPIDUtils = new AliEMCALPIDUtils();
91 
92  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),
141 {
142  for (Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ;
143  fMisalTransShift[i] = reco.fMisalTransShift[i] ; }
144  for (Int_t i = 0; i < 10 ; i++){ fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
145  for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
146  for (Int_t j = 0; j < 5 ; j++) { fMCGenerToAccept[j] = reco.fMCGenerToAccept[j] ; }
147  for (Int_t j = 0; j < 4 ; j++) { fBadStatusSelection[j] = reco.fBadStatusSelection[j] ; }
148 
149  if(reco.fEMCALBadChannelMap) {
150  // Copy constructor - not taking ownership over calibration histograms
151  fEMCALBadChannelMap = new TObjArray(reco.fEMCALBadChannelMap->GetEntries());
152  fEMCALBadChannelMap->SetOwner(false);
153  for(int ism = 0; ism < reco.fEMCALBadChannelMap->GetEntries(); ism++) fEMCALBadChannelMap->AddAt(reco.fEMCALBadChannelMap->At(ism), ism);
154  }
155 
156  if(reco.fEMCALRecalibrationFactors) {
157  // Copy constructor - not taking ownership over calibration histograms
159  fEMCALRecalibrationFactors->SetOwner(false);
160  for(int ism = 0; ism < reco.fEMCALRecalibrationFactors->GetEntries(); ism++) fEMCALRecalibrationFactors->AddAt(reco.fEMCALRecalibrationFactors->At(ism), ism);
161  }
162 
164  // Copy constructor - not taking ownership over calibration histograms
166  fEMCALTimeRecalibrationFactors->SetOwner(false);
167  for(int ism = 0; ism < reco.fEMCALTimeRecalibrationFactors->GetEntries(); ism++) fEMCALTimeRecalibrationFactors->AddAt(reco.fEMCALTimeRecalibrationFactors->At(ism), ism);
168  }
169 }
170 
174 //______________________________________________________________________
176 {
177  if (this == &reco)return *this;
178  ((TNamed *)this)->operator=(reco);
179 
180  for (Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ;
181  fMisalRotShift[i] = reco.fMisalRotShift[i] ; }
182  for (Int_t i = 0; i < 10 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
183  for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
184 
186  fPosAlgo = reco.fPosAlgo;
187  fW0 = reco.fW0;
189 
193 
196 
199  fLowGain = reco.fLowGain;
200 
203 
205 
208 
211 
217 
218  fPIDUtils = reco.fPIDUtils;
219 
223 
226  fCutR = reco.fCutR;
227  fCutEta = reco.fCutEta;
228  fCutPhi = reco.fCutPhi;
230  fMass = reco.fMass;
231  fStepSurface = reco.fStepSurface;
232  fStepCluster = reco.fStepCluster;
233  fITSTrackSA = reco.fITSTrackSA;
235 
250 
253  for (Int_t j = 0; j < 5 ; j++)
254  fMCGenerToAccept[j] = reco.fMCGenerToAccept[j];
255 
256  //
257  // Assign or copy construct the different TArrays
258  //
259  if (reco.fResidualEta)
260  {
261  if (fResidualEta)
262  *fResidualEta = *reco.fResidualEta;
263  else
264  fResidualEta = new TArrayF(*reco.fResidualEta);
265  }
266  else
267  {
268  if (fResidualEta) delete fResidualEta;
269  fResidualEta = 0;
270  }
271 
272  if (reco.fResidualPhi)
273  {
274  if (fResidualPhi)
275  *fResidualPhi = *reco.fResidualPhi;
276  else
277  fResidualPhi = new TArrayF(*reco.fResidualPhi);
278  }
279  else
280  {
281  if (fResidualPhi) delete fResidualPhi;
282  fResidualPhi = 0;
283  }
284 
285  if (reco.fMatchedTrackIndex)
286  {
287  if (fMatchedTrackIndex)
289  else
291  }
292  else
293  {
295  fMatchedTrackIndex = 0;
296  }
297 
298  if (reco.fMatchedClusterIndex)
299  {
302  else
304  }
305  else
306  {
309  }
310 
311  for (Int_t j = 0; j < 4 ; j++)
313 
315  if(reco.fEMCALBadChannelMap) {
316  // Copy constructor - not taking ownership over calibration histograms
317  fEMCALBadChannelMap = new TObjArray(reco.fEMCALBadChannelMap->GetEntries());
318  fEMCALBadChannelMap->SetOwner(false);
319  for(int ism = 0; ism < reco.fEMCALBadChannelMap->GetEntries(); ism++) fEMCALBadChannelMap->AddAt(reco.fEMCALBadChannelMap->At(ism), ism);
320  }
321 
323  if(reco.fEMCALRecalibrationFactors) {
324  // Copy constructor - not taking ownership over calibration histograms
326  fEMCALRecalibrationFactors->SetOwner(false);
327  for(int ism = 0; ism < reco.fEMCALRecalibrationFactors->GetEntries(); ism++) fEMCALRecalibrationFactors->AddAt(reco.fEMCALRecalibrationFactors->At(ism), ism);
328  }
329 
332  // Copy constructor - not taking ownership over calibration histograms
334  fEMCALTimeRecalibrationFactors->SetOwner(false);
335  for(int ism = 0; ism < reco.fEMCALTimeRecalibrationFactors->GetEntries(); ism++) fEMCALTimeRecalibrationFactors->AddAt(reco.fEMCALTimeRecalibrationFactors->At(ism), ism);
336  }
337 
340  // Copy constructor - not taking ownership over calibration histograms
342  fEMCALL1PhaseInTimeRecalibration->SetOwner(false);
343  for(int ism = 0; ism < reco.fEMCALL1PhaseInTimeRecalibration->GetEntries(); ism++) fEMCALL1PhaseInTimeRecalibration->AddAt(reco.fEMCALL1PhaseInTimeRecalibration->At(ism), ism);
344  }
345 
346  return *this;
347 }
348 
352 //_____________________________________
354 {
356  {
358  }
359 
361  {
363  }
364 
366  {
368  }
369 
370  if (fEMCALBadChannelMap)
371  {
372  delete fEMCALBadChannelMap;
373  }
374 
375  delete fMatchedTrackIndex ;
376  delete fMatchedClusterIndex ;
377  delete fResidualEta ;
378  delete fResidualPhi ;
379  delete fPIDUtils ;
380 
381  InitTrackCuts();
382 }
383 
396 //_______________________________________________________________________________
398  Float_t & amp, Double_t & time,
399  AliVCaloCells* cells)
400 {
401  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
402 
403  if(!geom)
404  {
405  AliError("No instance of the geometry is available");
406  return kFALSE;
407  }
408 
409  if ( absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules() )
410  return kFALSE;
411 
412  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1, status=0;
413 
414  if (!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
415  {
416  // cell absID does not exist
417  amp=0; time = 1.e9;
418  return kFALSE;
419  }
420 
421  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
422 
423  // Do not include bad channels found in analysis,
425  {
426  Bool_t bad = GetEMCALChannelStatus(imod, ieta, iphi,status);
427 
428  if ( status > 0 )
429  AliDebug(1,Form("Channel absId %d, status %d, set as bad %d",absID, status, bad));
430 
431  if ( bad ) return kFALSE;
432  }
433 
434  //Recalibrate energy
435  amp = cells->GetCellAmplitude(absID);
437  amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
438 
439  // Recalibrate time
440  time = cells->GetCellTime(absID);
441  time-=fConstantTimeShift*1e-9; // only in case of old Run1 simulation
442  Bool_t isLowGain = !(cells->GetCellHighGain(absID));//HG = false -> LG = true
443 
444  RecalibrateCellTime(absID,bc,time,isLowGain);
445 
446  //Recalibrate time with L1 phase
447  RecalibrateCellTimeL1Phase(imod, bc, time);
448 
449  return kTRUE;
450 }
451 
461 //_____________________________________________________________________________
463  const AliVCluster* cluster,
464  AliVCaloCells* cells)
465 {
466  if (!cluster)
467  {
468  AliInfo("Cluster pointer null!");
469  return kFALSE;
470  }
471 
472  // If the distance to the border is 0 or negative just exit accept all clusters
473  if (cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 )
474  return kTRUE;
475 
476  Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
477  Bool_t shared = kFALSE;
478  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
479 
480  AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
481  absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
482 
483  if (absIdMax==-1) return kFALSE;
484 
485  // Check if the cell is close to the borders:
486  Bool_t okrow = kFALSE;
487  Bool_t okcol = kFALSE;
488 
489  if (iSM < 0 || iphi < 0 || ieta < 0 )
490  {
491  AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
492  iSM,ieta,iphi));
493  return kFALSE; // trick coverity
494  }
495 
496  // Check rows/phi
497  Int_t iPhiLast = 24;
498  if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_Half ) iPhiLast /= 2;
499  else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_3rd ) iPhiLast /= 3;// 1/3 sm case
500  else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Ext ) iPhiLast /= 3;// 1/3 sm case
501 
502  if(iphi >= fNCellsFromEMCALBorder && iphi < iPhiLast - fNCellsFromEMCALBorder) okrow = kTRUE;
503 
504  // Check columns/eta
505  Int_t iEtaLast = 48;
506  if ( !fNoEMCALBorderAtEta0 || geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Standard )
507  {
508  // consider inner border
509  if ( geom->IsDCALSM(iSM) ) iEtaLast = iEtaLast*2/3;
510 
511  if ( ieta > fNCellsFromEMCALBorder && ieta < iEtaLast-fNCellsFromEMCALBorder ) okcol = kTRUE;
512  }
513  else
514  {
515  if (iSM%2==0)
516  {
517  if (ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
518  }
519  else
520  {
521  if (ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE;
522  }
523  }//eta 0 not checked
524 
525  AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
526  fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
527 
528  if (okcol && okrow)
529  {
530  return kTRUE;
531  }
532  else
533  {
534  AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
535 
536  return kFALSE;
537  }
538 }
539 
550 //_______________________________________________________________________________
552  const UShort_t* cellList,
553  Int_t nCells)
554 {
555  if (!fRemoveBadChannels) return kFALSE;
556  if (!fEMCALBadChannelMap) return kFALSE;
557 
558  Int_t icol = -1;
559  Int_t irow = -1;
560  Int_t imod = -1;
561  for (Int_t iCell = 0; iCell<nCells; iCell++)
562  {
563  //Get the column and row
564  Int_t iTower = -1, iIphi = -1, iIeta = -1;
565  geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
566 
567  if (fEMCALBadChannelMap->GetEntries() <= imod) continue;
568 
569  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
570 
571  Int_t status = 0;
572  if (GetEMCALChannelStatus(imod, icol, irow, status))
573  {
574  AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d, status %d\n",imod, icol, irow, status));
575  return kTRUE;
576  }
577  }// cell cluster loop
578 
579  return kFALSE;
580 }
581 
593 //___________________________________________________________________________
595  AliVCaloCells* cells, Int_t bc)
596 {
597  AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
598 
599  if(!geom)
600  {
601  AliError("No instance of the geometry is available");
602  return -1;
603  }
604 
605  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
606  geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
607  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
608 
609  // Get close cells index, energy and time, not in corners
610 
611  Int_t absID1 = -1;
612  Int_t absID2 = -1;
613 
614  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
615  if ( iphi > 0 ) absID2 = geom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
616 
617  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
618 
619  Int_t absID3 = -1;
620  Int_t absID4 = -1;
621 
622  if ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) )
623  {
624  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
625  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
626  }
627  else if ( ieta == 0 && imod%2 )
628  {
629  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
630  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
631  }
632  else
633  {
634  if ( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
635  absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
636  if ( ieta > 0 )
637  absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
638  }
639 
640  //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4);
641 
642  Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
643  Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
644 
645  AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
646  AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
647  AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
648  AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
649 
650  if (TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
651  if (TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
652  if (TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
653  if (TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
654 
655  return ecell1+ecell2+ecell3+ecell4;
656 }
657 
668 //_____________________________________________________________________________________________
669 Bool_t AliEMCALRecoUtils::IsExoticCell(Int_t absID, AliVCaloCells* cells, Int_t bc)
670 {
671  if (!fRejectExoticCells) return kFALSE;
672 
673  Float_t ecell = 0;
674  Double_t tcell = 0;
675  Bool_t accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
676 
677  if (!accept) return kTRUE; // reject this cell
678 
679  if (ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
680 
681  Float_t eCross = GetECross(absID,tcell,cells,bc);
682 
683  if (1-eCross/ecell > fExoticCellFraction)
684  {
685  AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",
686  absID,ecell,eCross,1-eCross/ecell));
687  return kTRUE;
688  }
689 
690  return kFALSE;
691 }
692 
702 //___________________________________________________________________
703 Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster,
704  AliVCaloCells *cells,
705  Int_t bc)
706 {
707  if (!cluster)
708  {
709  AliInfo("Cluster pointer null!");
710  return kFALSE;
711  }
712 
713  if (!fRejectExoticCluster) return kFALSE;
714 
715  // Get highest energy tower
716  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
717 
718  if(!geom)
719  {
720  AliError("No instance of the geometry is available");
721  return kFALSE;
722  }
723 
724  Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
725  Bool_t shared = kFALSE;
726  GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
727 
728  return IsExoticCell(absId,cells,bc);
729 }
730 
739 //_______________________________________________________________________
740 Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster)
741 {
742  if (!cluster)
743  {
744  AliInfo("Cluster pointer null!");
745  return 0;
746  }
747 
748  Float_t energy = cluster->E() ;
749  Float_t rdmEnergy = energy ;
750 
751  if (fSmearClusterEnergy)
752  {
753  rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
754  fSmearClusterParam[1] * energy +
755  fSmearClusterParam[2] );
756  AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
757  }
758 
759  return rdmEnergy;
760 }
761 
770 //____________________________________________________________________________
772 {
773  if (!cluster)
774  {
775  AliInfo("Cluster pointer null!");
776  return 0;
777  }
778 
779  Float_t energy = cluster->E();
780  if (energy < 0.100)
781  {
782  // Clusters with less than 100 MeV or negative are not possible
783  AliInfo(Form("Too low cluster energy!, E = %f < 0.100 GeV",energy));
784  return 0;
785  }
786  else if(energy > 300.)
787  {
788  AliInfo(Form("Too high cluster energy!, E = %f GeV, do not apply non linearity",energy));
789  return energy;
790  }
791 
792  switch (fNonLinearityFunction)
793  {
794  case kPi0MC:
795  {
796  //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
797  //fNonLinearityParams[0] = 1.014;
798  //fNonLinearityParams[1] =-0.03329;
799  //fNonLinearityParams[2] =-0.3853;
800  //fNonLinearityParams[3] = 0.5423;
801  //fNonLinearityParams[4] =-0.4335;
802  energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
803  ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
804  exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
805  break;
806  }
807 
808  case kPi0MCv2:
809  {
810  //Non-Linearity correction (from MC with function [0]/((x+[1])^[2]))+1;
811  //fNonLinearityParams[0] = 3.11111e-02;
812  //fNonLinearityParams[1] =-5.71666e-02;
813  //fNonLinearityParams[2] = 5.67995e-01;
814 
815  energy *= fNonLinearityParams[0]/TMath::Power(energy+fNonLinearityParams[1],fNonLinearityParams[2])+1;
816  break;
817  }
818 
819  case kPi0MCv3:
820  {
821  //Same as beam test corrected, change parameters
822  //fNonLinearityParams[0] = 9.81039e-01
823  //fNonLinearityParams[1] = 1.13508e-01;
824  //fNonLinearityParams[2] = 1.00173e+00;
825  //fNonLinearityParams[3] = 9.67998e-02;
826  //fNonLinearityParams[4] = 2.19381e+02;
827  //fNonLinearityParams[5] = 6.31604e+01;
828  //fNonLinearityParams[6] = 1;
829  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
830 
831  break;
832  }
833 
834 
835  case kPi0GammaGamma:
836  {
837  //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
838  //fNonLinearityParams[0] = 1.04;
839  //fNonLinearityParams[1] = -0.1445;
840  //fNonLinearityParams[2] = 1.046;
841  energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
842  break;
843  }
844 
845  case kPi0GammaConversion:
846  {
847  //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
848  //fNonLinearityParams[0] = 0.139393/0.1349766;
849  //fNonLinearityParams[1] = 0.0566186;
850  //fNonLinearityParams[2] = 0.982133;
851  energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
852 
853  break;
854  }
855 
856  case kBeamTest:
857  {
858  //From beam test, Alexei's results, for different ZS thresholds
859  // th=30 MeV; th = 45 MeV; th = 75 MeV
860  //fNonLinearityParams[0] = 1.007; 1.003; 1.002
861  //fNonLinearityParams[1] = 0.894; 0.719; 0.797
862  //fNonLinearityParams[2] = 0.246; 0.334; 0.358
863  //Rescale the param[0] with 1.03
864  energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
865 
866  break;
867  }
868 
869  case kBeamTestCorrected:
870  {
871  // From beam test, corrected for material between beam and EMCAL
872  //fNonLinearityParams[0] = 0.99078
873  //fNonLinearityParams[1] = 0.161499;
874  //fNonLinearityParams[2] = 0.655166;
875  //fNonLinearityParams[3] = 0.134101;
876  //fNonLinearityParams[4] = 163.282;
877  //fNonLinearityParams[5] = 23.6904;
878  //fNonLinearityParams[6] = 0.978;
879  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
880 
881  break;
882  }
883 
885  {
886  // From beam test, corrected for material between beam and EMCAL
887  // Different function to kBeamTestCorrected
888  //fNonLinearityParams[0] = 0.983504;
889  //fNonLinearityParams[1] = 0.210106;
890  //fNonLinearityParams[2] = 0.897274;
891  //fNonLinearityParams[3] = 0.0829064;
892  //fNonLinearityParams[4] = 152.299;
893  //fNonLinearityParams[5] = 31.5028;
894  //fNonLinearityParams[6] = 0.968;
895  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
896 
897  break;
898  }
899 
901  {
902  // Same function as kBeamTestCorrectedv2, different default parametrization.
903  //fNonLinearityParams[0] = 0.976941;
904  //fNonLinearityParams[1] = 0.162310;
905  //fNonLinearityParams[2] = 1.08689;
906  //fNonLinearityParams[3] = 0.0819592;
907  //fNonLinearityParams[4] = 152.338;
908  //fNonLinearityParams[5] = 30.9594;
909  //fNonLinearityParams[6] = 0.9615;
910  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
911 
912  break;
913  }
914 
915  case kSDMv5:
916  {
917  //Based on fit to the MC/data using kNoCorrection on the data - utilizes symmetric decay method and kPi0MCv5(MC) - 28 Oct 2013
918  //fNonLinearityParams[0] = 1.0;
919  //fNonLinearityParams[1] = 6.64778e-02;
920  //fNonLinearityParams[2] = 1.570;
921  //fNonLinearityParams[3] = 9.67998e-02;
922  //fNonLinearityParams[4] = 2.19381e+02;
923  //fNonLinearityParams[5] = 6.31604e+01;
924  //fNonLinearityParams[6] = 1.01286;
925  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));
926 
927  break;
928  }
929 
930  case kPi0MCv5:
931  {
932  //Based on comparing MC truth information to the reconstructed energy of clusters.
933  //fNonLinearityParams[0] = 1.0;
934  //fNonLinearityParams[1] = 6.64778e-02;
935  //fNonLinearityParams[2] = 1.570;
936  //fNonLinearityParams[3] = 9.67998e-02;
937  //fNonLinearityParams[4] = 2.19381e+02;
938  //fNonLinearityParams[5] = 6.31604e+01;
939  //fNonLinearityParams[6] = 1.01286;
940  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
941 
942  break;
943  }
944 
945  case kSDMv6:
946  {
947  //Based on fit to the MC/data using kNoCorrection on the data
948  // - utilizes symmetric decay method and kPi0MCv6(MC) - 09 Dec 2014
949  // - parameters constrained by the test beam data as well
950  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
951  //fNonLinearityParams[0] = 1.0;
952  //fNonLinearityParams[1] = 0.237767;
953  //fNonLinearityParams[2] = 0.651203;
954  //fNonLinearityParams[3] = 0.183741;
955  //fNonLinearityParams[4] = 155.427;
956  //fNonLinearityParams[5] = 17.0335;
957  //fNonLinearityParams[6] = 0.987054;
958  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
959 
960  break;
961  }
962 
963  case kPi0MCv6:
964  {
965  //Based on comparing MC truth information to the reconstructed energy of clusters.
966  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
967  //fNonLinearityParams[0] = 1.0;
968  //fNonLinearityParams[1] = 0.0797873;
969  //fNonLinearityParams[2] = 1.68322;
970  //fNonLinearityParams[3] = 0.0806098;
971  //fNonLinearityParams[4] = 244.586;
972  //fNonLinearityParams[5] = 116.938;
973  //fNonLinearityParams[6] = 1.00437;
974  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
975 
976  break;
977  }
978 
979  case kPCMv1:
980  {
981  //based on symmetric decays of pi0 meson
982  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
983  // parameters vary from MC to MC
984  //fNonLinearityParams[0] = 0.984876;
985  //fNonLinearityParams[1] = -9.999609;
986  //fNonLinearityParams[2] = -4.999891;
987  //fNonLinearityParams[3] = 0.;
988  //fNonLinearityParams[4] = 0.;
989  //fNonLinearityParams[5] = 0.;
990  //fNonLinearityParams[6] = 0.;
991  energy /= TMath::Power(fNonLinearityParams[0] + exp(fNonLinearityParams[1] + fNonLinearityParams[2]*energy),2);//result coming from calo-conv method
992 
993  break;
994  }
995 
996  case kPCMplusBTCv1:
997  {
998  //convolution of TestBeamCorrectedv3 with PCM method
999  //Based on comparing MC truth information to the reconstructed energy of clusters.
1000  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1001  // parameters vary from MC to MC
1002  //fNonLinearityParams[0] = 0.976941;
1003  //fNonLinearityParams[1] = 0.162310;
1004  //fNonLinearityParams[2] = 1.08689;
1005  //fNonLinearityParams[3] = 0.0819592;
1006  //fNonLinearityParams[4] = 152.338;
1007  //fNonLinearityParams[5] = 30.9594;
1008  //fNonLinearityParams[6] = 0.9615;
1009  //fNonLinearityParams[7] = 0.984876;
1010  //fNonLinearityParams[8] = -9.999609;
1011  //fNonLinearityParams[9] = -4.999891;
1012  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1013  energy /= TMath::Power(fNonLinearityParams[7] + exp(fNonLinearityParams[8] + fNonLinearityParams[9]*energy),2);//result coming from calo-conv method
1014 
1015  break;
1016  }
1017 
1018  case kPCMsysv1:
1019  {
1020  // Systematic variation of kPCMv1
1021  //Based on comparing MC truth information to the reconstructed energy of clusters.
1022  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1023  // parameters vary from MC to MC
1024  //fNonLinearityParams[0] = 0.0;
1025  //fNonLinearityParams[1] = 1.0;
1026  //fNonLinearityParams[2] = 1.0;
1027  //fNonLinearityParams[3] = 0.0;
1028  //fNonLinearityParams[4] = 1.0;
1029  //fNonLinearityParams[5] = 0.0;
1030  //fNonLinearityParams[6] = 0.0;
1031  energy /= TMath::Power( (fNonLinearityParams[0] + fNonLinearityParams[1] * TMath::Power(energy,fNonLinearityParams[2]) ) /
1032  (fNonLinearityParams[3] + fNonLinearityParams[4] * TMath::Power(energy,fNonLinearityParams[5]) ) + fNonLinearityParams[6], 2);//result coming from calo-conv method
1033 
1034  break;
1035  }
1036 
1037 
1038  case kNoCorrection:
1039  AliDebug(2,"No correction on the energy\n");
1040  break;
1041 
1042  }
1043 
1044  return energy;
1045 }
1046 
1051 //__________________________________________________
1053 {
1054  if (fNonLinearityFunction == kPi0MC) {
1055  fNonLinearityParams[0] = 1.014;
1056  fNonLinearityParams[1] = -0.03329;
1057  fNonLinearityParams[2] = -0.3853;
1058  fNonLinearityParams[3] = 0.5423;
1059  fNonLinearityParams[4] = -0.4335;
1060  }
1061 
1063  fNonLinearityParams[0] = 3.11111e-02;
1064  fNonLinearityParams[1] =-5.71666e-02;
1065  fNonLinearityParams[2] = 5.67995e-01;
1066  }
1067 
1069  fNonLinearityParams[0] = 9.81039e-01;
1070  fNonLinearityParams[1] = 1.13508e-01;
1071  fNonLinearityParams[2] = 1.00173e+00;
1072  fNonLinearityParams[3] = 9.67998e-02;
1073  fNonLinearityParams[4] = 2.19381e+02;
1074  fNonLinearityParams[5] = 6.31604e+01;
1075  fNonLinearityParams[6] = 1;
1076  }
1077 
1079  fNonLinearityParams[0] = 1.04;
1080  fNonLinearityParams[1] = -0.1445;
1081  fNonLinearityParams[2] = 1.046;
1082  }
1083 
1085  fNonLinearityParams[0] = 0.139393;
1086  fNonLinearityParams[1] = 0.0566186;
1087  fNonLinearityParams[2] = 0.982133;
1088  }
1089 
1091  if (fNonLinearThreshold == 30) {
1092  fNonLinearityParams[0] = 1.007;
1093  fNonLinearityParams[1] = 0.894;
1094  fNonLinearityParams[2] = 0.246;
1095  }
1096  if (fNonLinearThreshold == 45) {
1097  fNonLinearityParams[0] = 1.003;
1098  fNonLinearityParams[1] = 0.719;
1099  fNonLinearityParams[2] = 0.334;
1100  }
1101  if (fNonLinearThreshold == 75) {
1102  fNonLinearityParams[0] = 1.002;
1103  fNonLinearityParams[1] = 0.797;
1104  fNonLinearityParams[2] = 0.358;
1105  }
1106  }
1107 
1109  fNonLinearityParams[0] = 0.99078;
1110  fNonLinearityParams[1] = 0.161499;
1111  fNonLinearityParams[2] = 0.655166;
1112  fNonLinearityParams[3] = 0.134101;
1113  fNonLinearityParams[4] = 163.282;
1114  fNonLinearityParams[5] = 23.6904;
1115  fNonLinearityParams[6] = 0.978;
1116  }
1117 
1119  // Parameters until November 2015, use now kBeamTestCorrectedv3
1120  fNonLinearityParams[0] = 0.983504;
1121  fNonLinearityParams[1] = 0.210106;
1122  fNonLinearityParams[2] = 0.897274;
1123  fNonLinearityParams[3] = 0.0829064;
1124  fNonLinearityParams[4] = 152.299;
1125  fNonLinearityParams[5] = 31.5028;
1126  fNonLinearityParams[6] = 0.968;
1127  }
1128 
1130 
1131  // New parametrization of kBeamTestCorrectedv2
1132  // excluding point at 0.5 GeV from Beam Test Data
1133  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1134 
1135  fNonLinearityParams[0] = 0.976941;
1136  fNonLinearityParams[1] = 0.162310;
1137  fNonLinearityParams[2] = 1.08689;
1138  fNonLinearityParams[3] = 0.0819592;
1139  fNonLinearityParams[4] = 152.338;
1140  fNonLinearityParams[5] = 30.9594;
1141  fNonLinearityParams[6] = 0.9615;
1142  }
1143 
1144  if (fNonLinearityFunction == kSDMv5) {
1145  fNonLinearityParams[0] = 1.0;
1146  fNonLinearityParams[1] = 6.64778e-02;
1147  fNonLinearityParams[2] = 1.570;
1148  fNonLinearityParams[3] = 9.67998e-02;
1149  fNonLinearityParams[4] = 2.19381e+02;
1150  fNonLinearityParams[5] = 6.31604e+01;
1151  fNonLinearityParams[6] = 1.01286;
1152  }
1153 
1155  fNonLinearityParams[0] = 1.0;
1156  fNonLinearityParams[1] = 6.64778e-02;
1157  fNonLinearityParams[2] = 1.570;
1158  fNonLinearityParams[3] = 9.67998e-02;
1159  fNonLinearityParams[4] = 2.19381e+02;
1160  fNonLinearityParams[5] = 6.31604e+01;
1161  fNonLinearityParams[6] = 1.01286;
1162  }
1163 
1164  if (fNonLinearityFunction == kSDMv6) {
1165  fNonLinearityParams[0] = 1.0;
1166  fNonLinearityParams[1] = 0.237767;
1167  fNonLinearityParams[2] = 0.651203;
1168  fNonLinearityParams[3] = 0.183741;
1169  fNonLinearityParams[4] = 155.427;
1170  fNonLinearityParams[5] = 17.0335;
1171  fNonLinearityParams[6] = 0.987054;
1172  }
1173 
1175  fNonLinearityParams[0] = 1.0;
1176  fNonLinearityParams[1] = 0.0797873;
1177  fNonLinearityParams[2] = 1.68322;
1178  fNonLinearityParams[3] = 0.0806098;
1179  fNonLinearityParams[4] = 244.586;
1180  fNonLinearityParams[5] = 116.938;
1181  fNonLinearityParams[6] = 1.00437;
1182  }
1183 
1184 if (fNonLinearityFunction == kPCMv1) {
1185  //parameters change from MC production to MC production, they need to set for each period
1186  fNonLinearityParams[0] = 0.984876;
1187  fNonLinearityParams[1] = -9.999609;
1188  fNonLinearityParams[2] = -4.999891;
1189  fNonLinearityParams[3] = 0.;
1190  fNonLinearityParams[4] = 0.;
1191  fNonLinearityParams[5] = 0.;
1192  fNonLinearityParams[6] = 0.;
1193  }
1194 
1196  // test beam corrected values convoluted with symmetric meson decays values
1197  // for test beam:
1198  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1199  // for PCM method:
1200  // https://aliceinfo.cern.ch/Notes/node/211
1201  fNonLinearityParams[0] = 0.976941;
1202  fNonLinearityParams[1] = 0.162310;
1203  fNonLinearityParams[2] = 1.08689;
1204  fNonLinearityParams[3] = 0.0819592;
1205  fNonLinearityParams[4] = 152.338;
1206  fNonLinearityParams[5] = 30.9594;
1207  fNonLinearityParams[6] = 0.9615;
1208  fNonLinearityParams[7] = 0.984876;
1209  fNonLinearityParams[8] = -9.999609;
1210  fNonLinearityParams[9] = -4.999891;
1211  }
1212 
1214  //systematics for kPCMv1
1215  // for PCM method:
1216  // https://aliceinfo.cern.ch/Notes/node/211
1217  fNonLinearityParams[0] = 0.0;
1218  fNonLinearityParams[1] = 1.0;
1219  fNonLinearityParams[2] = 1.0;
1220  fNonLinearityParams[3] = 0.0;
1221  fNonLinearityParams[4] = 1.0;
1222  fNonLinearityParams[5] = 0.0;
1223  fNonLinearityParams[6] = 0.0;
1224  }
1225 }
1226 
1234 //____________________________________________________________________
1236 {
1237  fBadStatusSelection[0] = all; // declare all as bad if true, never mind the other settings
1238  fBadStatusSelection[1] = dead;
1239  fBadStatusSelection[2] = hot;
1240  fBadStatusSelection[3] = warm;
1241 }
1242 
1253 //____________________________________________________________________
1255 {
1256  if(fEMCALBadChannelMap)
1257  status = (Int_t) ((TH2I*)fEMCALBadChannelMap->At(iSM))->GetBinContent(iCol,iRow);
1258  else
1259  status = 0; // Channel is ok by default
1260 
1261  if ( status == AliCaloCalibPedestal::kAlive )
1262  {
1263  return kFALSE; // Good channel
1264  }
1265  else
1266  {
1267  if ( fBadStatusSelection[0] == kTRUE )
1268  {
1269  return kTRUE; // consider bad hot, dead and warm
1270  }
1271  else
1272  {
1273  if ( fBadStatusSelection[AliCaloCalibPedestal::kDead] == kTRUE &&
1274  status == AliCaloCalibPedestal::kDead )
1275  return kTRUE; // consider bad dead
1276  else if ( fBadStatusSelection[AliCaloCalibPedestal::kHot] == kTRUE &&
1277  status == AliCaloCalibPedestal::kHot )
1278  return kTRUE; // consider bad hot
1279  else if ( fBadStatusSelection[AliCaloCalibPedestal::kWarning] == kTRUE &&
1280  status == AliCaloCalibPedestal::kWarning )
1281  return kTRUE; // consider bad warm
1282  }
1283  }
1284 
1285  AliWarning(Form("Careful, bad channel selection not properly done: ism %d, icol %d, irow %d, status %d,\n"
1286  " fBadAll %d, fBadHot %d, fBadWarm %d, fBadDead %d",
1287  iSM, iCol, iRow, status,
1290 
1291  return kFALSE; // if everything fails, accept it.
1292 }
1293 
1307 //____________________________________________________________________
1308 void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
1309  AliVCaloCells* cells,
1310  const AliVCluster* clu,
1311  Int_t & absId,
1312  Int_t & iSupMod,
1313  Int_t & ieta,
1314  Int_t & iphi,
1315  Bool_t & shared)
1316 {
1317  Double_t eMax = -1.;
1318  Double_t eCell = -1.;
1319  Float_t fraction = 1.;
1320  Float_t recalFactor = 1.;
1321  Int_t cellAbsId = -1 ;
1322 
1323  Int_t iTower = -1;
1324  Int_t iIphi = -1;
1325  Int_t iIeta = -1;
1326  Int_t iSupMod0= -1;
1327 
1328  if (!clu)
1329  {
1330  AliInfo("Cluster pointer null!");
1331  absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
1332  return;
1333  }
1334 
1335  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1336  {
1337  cellAbsId = clu->GetCellAbsId(iDig);
1338  fraction = clu->GetCellAmplitudeFraction(iDig);
1339  //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
1340  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1341 
1342  geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1343  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1344 
1345  if (iDig==0)
1346  {
1347  iSupMod0=iSupMod;
1348  } else if (iSupMod0!=iSupMod)
1349  {
1350  shared = kTRUE;
1351  //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
1352  }
1353 
1355  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1356 
1357  eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
1358  //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
1359  if (eCell > eMax)
1360  {
1361  eMax = eCell;
1362  absId = cellAbsId;
1363  //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
1364  }
1365  }// cell loop
1366 
1367  //Get from the absid the supermodule, tower and eta/phi numbers
1368  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1369  //Gives SuperModule and Tower numbers
1370  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1371  //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
1372  //printf("Max end---\n");
1373 }
1374 
1383 //_________________________________________________________
1385 {
1386  if (eCell > 0 && eCluster > 0)
1387  {
1388  if ( fW0 > 0 ) return TMath::Max( 0., fW0 + TMath::Log( eCell / eCluster ) ) ;
1389  else return TMath::Log( eCluster / eCell ) ;
1390  }
1391  else
1392  return 0. ;
1393 }
1394 
1398 //______________________________________
1400 {
1401  fParticleType = kPhoton;
1402  fPosAlgo = kUnchanged;
1403  fW0 = 4.5;
1404 
1406  fNonLinearThreshold = 30;
1407 
1408  fExoticCellFraction = 0.97;
1409  fExoticCellDiffTime = 1e6;
1411 
1412  fAODFilterMask = 128;
1413  fAODHybridTracks = kFALSE;
1414  fAODTPCOnlyTracks = kTRUE;
1415 
1416  fCutEtaPhiSum = kTRUE;
1417  fCutEtaPhiSeparate = kFALSE;
1418 
1419  fCutR = 0.05;
1420  fCutEta = 0.025;
1421  fCutPhi = 0.05;
1422 
1423  fClusterWindow = 100;
1424  fMass = 0.139;
1425 
1426  fStepSurface = 20.;
1427  fStepCluster = 5.;
1429 
1430  fCutMinTrackPt = 0;
1431  fCutMinNClusterTPC = -1;
1432  fCutMinNClusterITS = -1;
1433 
1434  fCutMaxChi2PerClusterTPC = 1e10;
1435  fCutMaxChi2PerClusterITS = 1e10;
1436 
1437  fCutRequireTPCRefit = kFALSE;
1438  fCutRequireITSRefit = kFALSE;
1439  fCutAcceptKinkDaughters = kFALSE;
1440 
1441  fCutMaxDCAToVertexXY = 1e10;
1442  fCutMaxDCAToVertexZ = 1e10;
1443  fCutDCAToVertex2D = kFALSE;
1444 
1445  fCutRequireITSStandAlone = kFALSE; //MARCEL
1446  fCutRequireITSpureSA = kFALSE; //Marcel
1447 
1448  // Misalignment matrices
1449  for (Int_t i = 0; i < 15 ; i++)
1450  {
1451  fMisalTransShift[i] = 0.;
1452  fMisalRotShift[i] = 0.;
1453  }
1454 
1455  // Non linearity
1456  for (Int_t i = 0; i < 10 ; i++) fNonLinearityParams[i] = 0.;
1457 
1458  // For kBeamTestCorrectedv2 case, but default is no correction
1459  fNonLinearityParams[0] = 0.983504;
1460  fNonLinearityParams[1] = 0.210106;
1461  fNonLinearityParams[2] = 0.897274;
1462  fNonLinearityParams[3] = 0.0829064;
1463  fNonLinearityParams[4] = 152.299;
1464  fNonLinearityParams[5] = 31.5028;
1465  fNonLinearityParams[6] = 0.968;
1466 
1467  // Cluster energy smearing
1468  fSmearClusterEnergy = kFALSE;
1469  fSmearClusterParam[0] = 0.07; // * sqrt E term
1470  fSmearClusterParam[1] = 0.00; // * E term
1471  fSmearClusterParam[2] = 0.00; // constant
1472 }
1473 
1477 //_____________________________________________________
1479 {
1480  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1481 
1482  // In order to avoid rewriting the same histograms
1483  Bool_t oldStatus = TH1::AddDirectoryStatus();
1484  TH1::AddDirectory(kFALSE);
1485 
1487  for (int i = 0; i < 22; i++)
1488  fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1489  Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
1490  //Init the histograms with 1
1491  for (Int_t sm = 0; sm < 22; sm++)
1492  {
1493  for (Int_t i = 0; i < 48; i++)
1494  {
1495  for (Int_t j = 0; j < 24; j++)
1496  {
1498  }
1499  }
1500  }
1501 
1502  fEMCALRecalibrationFactors->SetOwner(kTRUE);
1503  fEMCALRecalibrationFactors->Compress();
1504 
1505  // In order to avoid rewriting the same histograms
1506  TH1::AddDirectory(oldStatus);
1507 }
1508 
1512 //_________________________________________________________
1514 {
1515  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1516 
1517  // In order to avoid rewriting the same histograms
1518  Bool_t oldStatus = TH1::AddDirectoryStatus();
1519  TH1::AddDirectory(kFALSE);
1520 
1523 
1524  for (int i = 0; i < 4; i++)
1525  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1526  Form("hAllTimeAvBC%d",i),
1527  48*24*22,0.,48*24*22) );
1528  // Init the histograms with 0
1529  for (Int_t iBC = 0; iBC < 4; iBC++)
1530  {
1531  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1532  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kFALSE);
1533  }
1534 
1535  if(fLowGain) {
1536  for (int iBC = 0; iBC < 4; iBC++) {
1537  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvLGBC%d",iBC),
1538  Form("hAllTimeAvLGBC%d",iBC),
1539  48*24*22,0.,48*24*22) );
1540  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1541  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kTRUE);
1542  }
1543  }
1544 
1545 
1546 
1547  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1548  fEMCALTimeRecalibrationFactors->Compress();
1549 
1550  // In order to avoid rewriting the same histograms
1551  TH1::AddDirectory(oldStatus);
1552 }
1553 
1557 //____________________________________________________
1559 {
1560  AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1561 
1562  // In order to avoid rewriting the same histograms
1563  Bool_t oldStatus = TH1::AddDirectoryStatus();
1564  TH1::AddDirectory(kFALSE);
1565 
1566  fEMCALBadChannelMap = new TObjArray(22);
1567  //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
1568 
1569  for (int i = 0; i < 22; i++)
1570  fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1571 
1572  fEMCALBadChannelMap->SetOwner(kTRUE);
1573  fEMCALBadChannelMap->Compress();
1574 
1575  // In order to avoid rewriting the same histograms
1576  TH1::AddDirectory(oldStatus);
1577 }
1578 
1582 //___________________________________________________________
1584 {
1585  AliDebug(2,"AliEMCALRecoUtils::InitEMCALL1PhaseInTimeRecalibrationFactors()");
1586 
1587  // In order to avoid rewriting the same histograms
1588  Bool_t oldStatus = TH1::AddDirectoryStatus();
1589  TH1::AddDirectory(kFALSE);
1590 
1592 
1593  fEMCALL1PhaseInTimeRecalibration->Add(new TH1C("h0","EMCALL1phaseForSM", 22, 0, 22));
1594 
1595  for (Int_t i = 0; i < 22; i++) //loop over SMs, default value = 0
1597 
1598  fEMCALL1PhaseInTimeRecalibration->SetOwner(kTRUE);
1600 
1601  // In order to avoid rewriting the same histograms
1602  TH1::AddDirectory(oldStatus);
1603 }
1604 
1614 //____________________________________________________________________________
1615 void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
1616  AliVCluster * cluster,
1617  AliVCaloCells * cells,
1618  Int_t bc)
1619 {
1620  if (!cluster)
1621  {
1622  AliInfo("Cluster pointer null!");
1623  return;
1624  }
1625 
1626  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1627  UShort_t * index = cluster->GetCellsAbsId() ;
1628  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1629  Int_t ncells = cluster->GetNCells();
1630 
1631  // Initialize some used variables
1632  Float_t energy = 0;
1633  Int_t absId =-1;
1634  Int_t icol =-1, irow =-1, imod=1;
1635  Float_t factor = 1, frac = 0;
1636  Int_t absIdMax = -1;
1637  Float_t emax = 0;
1638 
1639  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1640  for (Int_t icell = 0; icell < ncells; icell++)
1641  {
1642  absId = index[icell];
1643  frac = fraction[icell];
1644  if (frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
1645 
1647  {
1648  // Energy
1649  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1650  geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1651  if (fEMCALRecalibrationFactors->GetEntries() <= imod)
1652  continue;
1653  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1654  factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1655 
1656  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1657  imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1658 
1659  }
1660 
1661  energy += cells->GetCellAmplitude(absId)*factor*frac;
1662 
1663  if (emax < cells->GetCellAmplitude(absId)*factor*frac)
1664  {
1665  emax = cells->GetCellAmplitude(absId)*factor*frac;
1666  absIdMax = absId;
1667  }
1668  }
1669 
1670  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
1671 
1672  cluster->SetE(energy);
1673 
1674  // Recalculate time of cluster
1675  Double_t timeorg = cluster->GetTOF();
1676  Bool_t isLowGain = !(cells->GetCellHighGain(absIdMax));//HG = false -> LG = true
1677 
1678  Double_t time = cells->GetCellTime(absIdMax);
1680  RecalibrateCellTime(absIdMax,bc,time,isLowGain);
1681  time-=fConstantTimeShift*1e-9; // only in case of Run1 old simulations
1682 
1683  // Recalibrate time with L1 phase
1685  RecalibrateCellTimeL1Phase(imod, bc, time);
1686 
1687  cluster->SetTOF(time);
1688 
1689  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
1690 }
1691 
1699 //_______________________________________________________________________
1700 void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, Int_t bc)
1701 {
1703  return;
1704 
1705  if (!cells)
1706  {
1707  AliInfo("Cells pointer null!");
1708  return;
1709  }
1710 
1711  Short_t absId =-1;
1712  Bool_t accept = kFALSE;
1713  Float_t ecell = 0;
1714  Double_t tcell = 0;
1715  Double_t ecellin = 0;
1716  Double_t tcellin = 0;
1717  Int_t mclabel = -1;
1718  Double_t efrac = 0;
1719 
1720  Int_t nEMcell = cells->GetNumberOfCells() ;
1721  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1722  {
1723  cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
1724 
1725  accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
1726  if (!accept)
1727  {
1728  ecell = 0;
1729  tcell = -1;
1730  }
1731 
1732  // Set new values
1733  cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
1734  }
1735 
1736  fCellsRecalibrated = kTRUE;
1737 }
1738 
1747 //_______________________________________________________________________________________________________
1748 void AliEMCALRecoUtils::RecalibrateCellTime(Int_t absId, Int_t bc, Double_t & celltime, Bool_t isLGon) const
1749 {
1750  if (!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0) {
1751  if(fLowGain)
1752  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,isLGon)*1.e-9;
1753  else
1754  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,kFALSE)*1.e-9;
1755  }
1756 }
1757 
1765 //_______________________________________________________________________________________________________
1767 {
1768  if (!fCellsRecalibrated && IsL1PhaseInTimeRecalibrationOn() && bc >= 0)
1769  {
1770  bc=bc%4;
1771 
1772  Float_t offsetPerSM=0.;
1773  Int_t l1PhaseShift = GetEMCALL1PhaseInTimeRecalibrationForSM(iSM);
1774  Int_t l1Phase=l1PhaseShift & 3; //bit operation
1775 
1776  if(bc >= l1Phase)
1777  offsetPerSM = (bc - l1Phase)*25;
1778  else
1779  offsetPerSM = (bc - l1Phase + 4)*25;
1780 
1781  Int_t l1shiftOffset=l1PhaseShift>>2; //bit operation
1782  l1shiftOffset*=25;
1783 
1784  celltime -= offsetPerSM*1.e-9;
1785  celltime -= l1shiftOffset*1.e-9;
1786  }
1787 }
1788 
1804 //______________________________________________________________________________
1805 void AliEMCALRecoUtils::RecalculateCellLabelsRemoveAddedGenerator( Int_t absID, AliVCluster* clus, AliMCEvent* mc,
1806  Float_t & amp, TArrayI & labeArr, TArrayF & eDepArr) const
1807 {
1808  TString genName ;
1809  Float_t eDepFrac[4];
1810 
1811  Float_t edepTotFrac = 1;
1812  Bool_t found = kFALSE;
1813  Float_t ampOrg = amp;
1814 
1815  //
1816  // Get the energy deposition fraction from cluster.
1817  //
1818  for(Int_t icluscell = 0; icluscell < clus->GetNCells(); icluscell++ )
1819  {
1820  if(absID == clus->GetCellAbsId(icluscell))
1821  {
1822  clus->GetCellMCEdepFractionArray(icluscell,eDepFrac);
1823 
1824  found = kTRUE;
1825 
1826  break;
1827  }
1828  }
1829 
1830  if ( !found )
1831  {
1832  AliWarning(Form("Cell abs ID %d NOT found in cluster",absID));
1833  return;
1834  }
1835 
1836  //
1837  // Check if there is a particle contribution from a given generator name.
1838  // If it is not one of the selected generators,
1839  // remove the constribution from the cell.
1840  //
1841  if ( mc && fNMCGenerToAccept > 0 )
1842  {
1843  //printf("Accept contribution from generator? \n");
1844  for(UInt_t imc = 0; imc < 4; imc++)
1845  {
1846  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc )
1847  {
1848  mc->GetCocktailGenerator(clus->GetLabelAt(imc),genName);
1849 
1850  Bool_t generOK = kFALSE;
1851  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
1852  {
1853  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
1854  }
1855 
1856  if ( !generOK )
1857  {
1858  amp-=ampOrg*eDepFrac[imc];
1859 
1860  edepTotFrac-=eDepFrac[imc];
1861 
1862  eDepFrac[imc] = 0;
1863  }
1864 
1865  } // eDep > 0
1866  } // 4 possible loop
1867 
1868  } // accept at least one kind of generator
1869 
1870  //
1871  // Add MC label and Edep to corresponding array (to be used later in digits)
1872  //
1873  Int_t nLabels = 0;
1874  for(UInt_t imc = 0; imc < 4; imc++)
1875  {
1876  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc && edepTotFrac > 0 )
1877  {
1878  nLabels++;
1879 
1880  labeArr.Set(nLabels);
1881  labeArr.AddAt(clus->GetLabelAt(imc), nLabels-1);
1882 
1883  eDepArr.Set(nLabels);
1884  eDepArr.AddAt( (eDepFrac[imc]/edepTotFrac) * amp, nLabels-1);
1885  // use as deposited energy a fraction of the simulated energy (smeared and with noise)
1886  } // edep > 0
1887 
1888  } // mc cell label loop
1889 
1890  //
1891  // If no label found, reject cell
1892  // It can happen to have this case (4 MC labels per cell is not enough for some cases)
1893  // better to remove. To be treated carefully.
1894  //
1895  if ( nLabels == 0 ) amp = 0;
1896 
1897 }
1898 
1906 //______________________________________________________________________________
1907 void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
1908  AliVCaloCells* cells,
1909  AliVCluster* clu)
1910 {
1911  if (!clu)
1912  {
1913  AliInfo("Cluster pointer null!");
1914  return;
1915  }
1916 
1918  else if (fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1919  else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
1920 }
1921 
1930 //_____________________________________________________________________________________________
1932  AliVCaloCells* cells,
1933  AliVCluster* clu)
1934 {
1935  Double_t eCell = 0.;
1936  Float_t fraction = 1.;
1937  Float_t recalFactor = 1.;
1938 
1939  Int_t absId = -1;
1940  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1941  Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
1942  Float_t weight = 0., totalWeight=0.;
1943  Float_t newPos[3] = {-1.,-1.,-1.};
1944  Double_t pLocal[3], pGlobal[3];
1945  Bool_t shared = kFALSE;
1946 
1947  Float_t clEnergy = clu->E(); //Energy already recalibrated previously
1948  if (clEnergy <= 0) return;
1949 
1950  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
1951  Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1952 
1953  //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1954 
1955  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1956  {
1957  absId = clu->GetCellAbsId(iDig);
1958  fraction = clu->GetCellAmplitudeFraction(iDig);
1959  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1960 
1961  if (!fCellsRecalibrated)
1962  {
1963  geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
1964  geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
1965  if (IsRecalibrationOn()) {
1966  recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1967  }
1968  }
1969 
1970  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1971 
1972  weight = GetCellWeight(eCell,clEnergy);
1973  totalWeight += weight;
1974 
1975  geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
1976  //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
1977  geom->GetGlobal(pLocal,pGlobal,iSupModMax);
1978  //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1979 
1980  for (int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
1981  }// cell loop
1982 
1983  if (totalWeight>0)
1984  {
1985  for (int i=0; i<3; i++ ) newPos[i] /= totalWeight;
1986  }
1987 
1988  //Float_t pos[]={0,0,0};
1989  //clu->GetPosition(pos);
1990  //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
1991  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1992 
1993  if (iSupModMax > 1) { //sector 1
1994  newPos[0] +=fMisalTransShift[3];//-=3.093;
1995  newPos[1] +=fMisalTransShift[4];//+=6.82;
1996  newPos[2] +=fMisalTransShift[5];//+=1.635;
1997  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
1998  } else { //sector 0
1999  newPos[0] +=fMisalTransShift[0];//+=1.134;
2000  newPos[1] +=fMisalTransShift[1];//+=8.2;
2001  newPos[2] +=fMisalTransShift[2];//+=1.197;
2002  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
2003  }
2004  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
2005 
2006  clu->SetPosition(newPos);
2007 }
2008 
2017 //____________________________________________________________________________________________
2019  AliVCaloCells* cells,
2020  AliVCluster* clu)
2021 {
2022  Double_t eCell = 1.;
2023  Float_t fraction = 1.;
2024  Float_t recalFactor = 1.;
2025 
2026  Int_t absId = -1;
2027  Int_t iTower = -1;
2028  Int_t iIphi = -1, iIeta = -1;
2029  Int_t iSupMod = -1, iSupModMax = -1;
2030  Int_t iphi = -1, ieta =-1;
2031  Bool_t shared = kFALSE;
2032 
2033  Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
2034 
2035  if (clEnergy <= 0)
2036  return;
2037  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
2038  Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
2039 
2040  Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
2041  Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
2042  Int_t startingSM = -1;
2043 
2044  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
2045  {
2046  absId = clu->GetCellAbsId(iDig);
2047  fraction = clu->GetCellAmplitudeFraction(iDig);
2048  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2049 
2050  if (iDig==0) startingSM = iSupMod;
2051  else if (iSupMod != startingSM) areInSameSM = kFALSE;
2052 
2053  eCell = cells->GetCellAmplitude(absId);
2054 
2055  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2056  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
2057 
2058  if (!fCellsRecalibrated)
2059  {
2060  if (IsRecalibrationOn()) {
2061  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2062  }
2063  }
2064 
2065  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2066 
2067  weight = GetCellWeight(eCell,clEnergy);
2068  if (weight < 0) weight = 0;
2069  totalWeight += weight;
2070  weightedCol += ieta*weight;
2071  weightedRow += iphi*weight;
2072 
2073  //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
2074  }// cell loop
2075 
2076  Float_t xyzNew[]={-1.,-1.,-1.};
2077  if (areInSameSM == kTRUE)
2078  {
2079  //printf("In Same SM\n");
2080  weightedCol = weightedCol/totalWeight;
2081  weightedRow = weightedRow/totalWeight;
2082  geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2083  }
2084  else
2085  {
2086  //printf("In Different SM\n");
2087  geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2088  }
2089 
2090  clu->SetPosition(xyzNew);
2091 }
2092 
2100 //___________________________________________________________________________________________
2102  AliVCaloCells* cells,
2103  AliVCluster * cluster)
2104 {
2105  if (!fRecalDistToBadChannels) return;
2106 
2107  if (!cluster)
2108  {
2109  AliInfo("Cluster pointer null!");
2110  return;
2111  }
2112 
2113  // Get channels map of the supermodule where the cluster is.
2114  Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
2115  Bool_t shared = kFALSE;
2116  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
2117  TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
2118 
2119  Int_t dRrow, dRcol;
2120  Float_t minDist = 10000.;
2121  Float_t dist = 0.;
2122 
2123  // Loop on tower status map
2124  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2125  {
2126  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2127  {
2128  // Check if tower is bad.
2129  if (hMap->GetBinContent(icol,irow)==0) continue;
2130  //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
2131  // iSupMod,icol, irow, icolM,irowM);
2132 
2133  dRrow=TMath::Abs(irowM-irow);
2134  dRcol=TMath::Abs(icolM-icol);
2135  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2136  if (dist < minDist)
2137  {
2138  //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
2139  minDist = dist;
2140  }
2141  }
2142  }
2143 
2144  // In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
2145  if (shared)
2146  {
2147  TH2D* hMap2 = 0;
2148  Int_t iSupMod2 = -1;
2149 
2150  // The only possible combinations are (0,1), (2,3) ... (8,9)
2151  if (iSupMod%2) iSupMod2 = iSupMod-1;
2152  else iSupMod2 = iSupMod+1;
2153  hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
2154 
2155  // Loop on tower status map of second super module
2156  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2157  {
2158  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2159  {
2160  // Check if tower is bad.
2161  if (hMap2->GetBinContent(icol,irow)==0) continue;
2162 
2163  //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",
2164  // iSupMod2,icol, irow,iSupMod,icolM,irowM);
2165  dRrow=TMath::Abs(irow-irowM);
2166 
2167  if (iSupMod%2)
2168  dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
2169  else
2170  dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
2171 
2172  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2173  if (dist < minDist) minDist = dist;
2174  }
2175  }
2176  }// shared cluster in 2 SuperModules
2177 
2178  AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
2179  cluster->SetDistanceToBadChannel(minDist);
2180 }
2181 
2188 //__________________________________________________________________
2189 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
2190 {
2191  if (!cluster)
2192  {
2193  AliInfo("Cluster pointer null!");
2194  return;
2195  }
2196 
2197  if (cluster->GetM02() != 0)
2198  fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
2199 
2200  Float_t pidlist[AliPID::kSPECIESCN+1];
2201  for (Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
2202 
2203  cluster->SetPID(pidlist);
2204 }
2205 
2226 //___________________________________________________________________________________________________________________
2228  AliVCaloCells* cells, AliVCluster * cluster,
2229  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2230  Float_t & enAfterCuts, Float_t & l0, Float_t & l1,
2231  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2232  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2233 {
2234  if (!cluster)
2235  {
2236  AliInfo("Cluster pointer null!");
2237  return;
2238  }
2239 
2240  Double_t eCell = 0.;
2241  Double_t tCell = 0.;
2242  Float_t fraction = 1.;
2243  Float_t recalFactor = 1.;
2244  Bool_t isLowGain = kFALSE;
2245 
2246  Int_t iSupMod = -1;
2247  Int_t iTower = -1;
2248  Int_t iIphi = -1;
2249  Int_t iIeta = -1;
2250  Int_t iphi = -1;
2251  Int_t ieta = -1;
2252  Double_t etai = -1.;
2253  Double_t phii = -1.;
2254 
2255  Int_t nstat = 0 ;
2256  Float_t wtot = 0.;
2257  Double_t w = 0.;
2258  Double_t etaMean = 0.;
2259  Double_t phiMean = 0.;
2260 
2261  Double_t pGlobal[3];
2262 
2263  // Loop on cells, calculate the cluster energy, in case a cut on cell energy is added,
2264  // or the non linearity correction was applied
2265  // and to check if the cluster is between 2 SM in eta
2266  Int_t iSM0 = -1;
2267  Bool_t shared = kFALSE;
2268  Float_t energy = 0;
2269 
2270  enAfterCuts = 0;
2271  l0 = 0; l1 = 0;
2272  disp = 0; dEta = 0; dPhi = 0;
2273  sEta = 0; sPhi = 0; sEtaPhi = 0;
2274 
2275  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2276  {
2277  // Get from the absid the supermodule, tower and eta/phi numbers
2278 
2279  Int_t absId = cluster->GetCellAbsId(iDigit);
2280 
2281  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2282  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2283 
2284  // Check if there are cells of different SM
2285  if (iDigit == 0 ) iSM0 = iSupMod;
2286  else if (iSupMod!= iSM0) shared = kTRUE;
2287 
2288  // Get the cell energy, if recalibration is on, apply factors
2289  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2290  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2291 
2292  if (IsRecalibrationOn())
2293  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2294 
2295  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2296  tCell = cells->GetCellTime (absId);
2297  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2298 
2299  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2300  tCell*=1e9;
2301  tCell-=fConstantTimeShift;
2302 
2303  if(eCell > 0.05) // at least a minimum 50 MeV cut
2304  energy += eCell;
2305 
2306  if(eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut)
2307  enAfterCuts += eCell;
2308  } // cell loop
2309 
2310 
2311  // Loop on cells to calculate weights and shower shape terms parameters
2312  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2313  {
2314  // Get from the absid the supermodule, tower and eta/phi numbers
2315  Int_t absId = cluster->GetCellAbsId(iDigit);
2316 
2317  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2318  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2319 
2320  //Get the cell energy, if recalibration is on, apply factors
2321  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2322  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2323 
2325  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2326 
2327  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2328  tCell = cells->GetCellTime (absId);
2329  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2330 
2331  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2332  tCell*=1e9;
2333  tCell-=fConstantTimeShift;
2334 
2335  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2336  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2337  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2338 
2339  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2340  {
2341  w = GetCellWeight(eCell, energy);
2342 
2343  // Cell index
2344  if ( fShowerShapeCellLocationType == 0 )
2345  {
2346  etai=(Double_t)ieta;
2347  phii=(Double_t)iphi;
2348  }
2349  // Cell angle location
2350  else if( fShowerShapeCellLocationType == 1 )
2351  {
2352  geom->EtaPhiFromIndex(absId, etai, phii);
2353  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2354  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2355  }
2356  else
2357  {
2358  geom->GetGlobal(absId,pGlobal);
2359 
2360  // Cell x-z location
2361  if( fShowerShapeCellLocationType == 2 )
2362  {
2363  etai = pGlobal[2];
2364  phii = pGlobal[0];
2365  }
2366  // Cell r-z location
2367  else
2368  {
2369  etai = pGlobal[2];
2370  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2371  }
2372  }
2373 
2374  if (w > 0.0)
2375  {
2376  wtot += w ;
2377  nstat++;
2378 
2379  // Shower shape
2380  sEta += w * etai * etai ;
2381  etaMean += w * etai ;
2382  sPhi += w * phii * phii ;
2383  phiMean += w * phii ;
2384  sEtaPhi += w * etai * phii ;
2385  }
2386  }
2387  else if(eCell > 0.05)
2388  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2389  } // cell loop
2390 
2391  //printf("sEta %f sPhi %f etaMean %f phiMean %f sEtaPhi %f wtot %f\n",sEta,sPhi,etaMean,phiMean,sEtaPhi, wtot);
2392 
2393  // Normalize to the weight
2394  if (wtot > 0)
2395  {
2396  etaMean /= wtot ;
2397  phiMean /= wtot ;
2398  }
2399  else
2400  AliDebug(2,Form("Wrong weight %f\n", wtot));
2401 
2402  // Loop on cells to calculate dispersion
2403  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2404  {
2405  // Get from the absid the supermodule, tower and eta/phi numbers
2406  Int_t absId = cluster->GetCellAbsId(iDigit);
2407 
2408  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2409  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2410 
2411  //Get the cell energy, if recalibration is on, apply factors
2412  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2413  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2414  if (IsRecalibrationOn())
2415  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2416 
2417  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2418  tCell = cells->GetCellTime (absId);
2419  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2420 
2421  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2422  tCell*=1e9;
2423  tCell-=fConstantTimeShift;
2424 
2425  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2426  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2427  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2428 
2429  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2430  {
2431  w = GetCellWeight(eCell,cluster->E());
2432 
2433  // Cell index
2434  if ( fShowerShapeCellLocationType == 0 )
2435  {
2436  etai=(Double_t)ieta;
2437  phii=(Double_t)iphi;
2438  }
2439  // Cell angle location
2440  else if( fShowerShapeCellLocationType == 1 )
2441  {
2442  geom->EtaPhiFromIndex(absId, etai, phii);
2443  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2444  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2445  }
2446  else
2447  {
2448  geom->GetGlobal(absId,pGlobal);
2449 
2450  // Cell x-z location
2451  if( fShowerShapeCellLocationType == 2 )
2452  {
2453  etai = pGlobal[2];
2454  phii = pGlobal[0];
2455  }
2456  // Cell r-z location
2457  else
2458  {
2459  etai = pGlobal[2];
2460  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2461  }
2462  }
2463 
2464  if (w > 0.0)
2465  {
2466  disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
2467  dEta += w * (etai-etaMean)*(etai-etaMean) ;
2468  dPhi += w * (phii-phiMean)*(phii-phiMean) ;
2469  }
2470  }
2471  else
2472  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2473  }// cell loop
2474 
2475  // Normalize to the weigth and set shower shape parameters
2476  if (wtot > 0 && nstat > 1)
2477  {
2478  disp /= wtot ;
2479  dEta /= wtot ;
2480  dPhi /= wtot ;
2481  sEta /= wtot ;
2482  sPhi /= wtot ;
2483  sEtaPhi /= wtot ;
2484 
2485  sEta -= etaMean * etaMean ;
2486  sPhi -= phiMean * phiMean ;
2487  sEtaPhi -= etaMean * phiMean ;
2488 
2489  l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2490  l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2491 
2492  //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);
2493 
2494  }
2495  else
2496  {
2497  l0 = 0. ;
2498  l1 = 0. ;
2499  dEta = 0. ; dPhi = 0. ; disp = 0. ;
2500  sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
2501  }
2502 }
2503 
2522 //___________________________________________________________________________________________________________________
2524  AliVCaloCells* cells, AliVCluster * cluster,
2525  Float_t & l0, Float_t & l1,
2526  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2527  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2528 {
2529  Float_t newEnergy = 0;
2530  Float_t cellEmin = 0.05; // 50 MeV
2531  Float_t cellTimeWindow = 1000; // open cut
2532  Int_t bc = 0;
2534  cellEmin, cellTimeWindow, bc,
2535  newEnergy, l0, l1, disp,
2536  dEta, dPhi, sEta, sPhi, sEtaPhi);
2537 }
2538 
2547 //____________________________________________________________________________________________
2549  AliVCaloCells* cells,
2550  AliVCluster * cluster)
2551 {
2552  Float_t l0 = 0., l1 = 0.;
2553  Float_t disp = 0., dEta = 0., dPhi = 0.;
2554  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2555 
2556  AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
2557  dEta, dPhi, sEta, sPhi, sEtaPhi);
2558 
2559  cluster->SetM02(l0);
2560  cluster->SetM20(l1);
2561  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2562 
2563 }
2564 
2577 //____________________________________________________________________________________________
2579  AliVCaloCells* cells, AliVCluster * cluster,
2580  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2581  Float_t & enAfterCuts)
2582 {
2583  Float_t l0 = 0., l1 = 0.;
2584  Float_t disp = 0., dEta = 0., dPhi = 0.;
2585  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2586 
2588  cellEcut, cellTimeCut, bc,
2589  enAfterCuts, l0, l1, disp,
2590  dEta, dPhi, sEta, sPhi, sEtaPhi);
2591 
2592  cluster->SetM02(l0);
2593  cluster->SetM20(l1);
2594  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2595 
2596 }
2597 
2611 //____________________________________________________________________________
2612 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
2613  TObjArray * clusterArr,
2614  const AliEMCALGeometry *geom,
2615  AliMCEvent * mc)
2616 {
2617  fMatchedTrackIndex ->Reset();
2618  fMatchedClusterIndex->Reset();
2619  fResidualPhi->Reset();
2620  fResidualEta->Reset();
2621 
2622  fMatchedTrackIndex ->Set(1000);
2623  fMatchedClusterIndex->Set(1000);
2624  fResidualPhi->Set(1000);
2625  fResidualEta->Set(1000);
2626 
2627  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
2628  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
2629 
2630  // Init the magnetic field if not already on
2631  if (!TGeoGlobalMagField::Instance()->GetField())
2632  {
2633  if (!event->InitMagneticField())
2634  {
2635  AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
2636  }
2637  } // Init mag field
2638 
2639  if (esdevent)
2640  {
2641  UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
2642  UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
2643  Bool_t desc1 = (mask1 >> 3) & 0x1;
2644  Bool_t desc2 = (mask2 >> 3) & 0x1;
2645  if (desc1==0 || desc2==0) {
2646  // AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
2647  // mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
2648  // mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
2649  fITSTrackSA=kTRUE;
2650  }
2651  }
2652 
2653  TObjArray *clusterArray = 0x0;
2654  if (!clusterArr)
2655  {
2656  clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
2657  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2658  {
2659  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2660  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells()))
2661  continue;
2662  clusterArray->AddAt(cluster,icl);
2663  }
2664  }
2665 
2666  Int_t matched=0;
2667  Double_t cv[21];
2668  TString genName;
2669  for (Int_t i=0; i<21;i++) cv[i]=0;
2670  for (Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
2671  {
2672  AliExternalTrackParam *trackParam = 0;
2673  Int_t mcLabel = -1;
2674  // If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
2675  AliESDtrack *esdTrack = 0;
2676  AliAODTrack *aodTrack = 0;
2677  if (esdevent)
2678  {
2679  esdTrack = esdevent->GetTrack(itr);
2680  if (!esdTrack) continue;
2681  if (!IsAccepted(esdTrack)) continue;
2682  if (esdTrack->Pt()<fCutMinTrackPt) continue;
2683 
2684  if ( TMath::Abs(esdTrack->Eta()) > 0.9 ) continue;
2685 
2686  // Save some time and memory in case of no DCal present
2687  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2688  {
2689  Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
2690  if ( phi <= 10 || phi >= 250 ) continue;
2691  }
2692 
2693  if (!fITSTrackSA)
2694  trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); // if TPC Available
2695  else
2696  trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
2697 
2698  mcLabel = TMath::Abs(esdTrack->GetLabel());
2699  }
2700 
2701  // If the input event is AOD, the starting point for extrapolation is at vertex
2702  // AOD tracks are selected according to its filterbit.
2703  else if (aodevent)
2704  {
2705  aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(itr));
2706 
2707  if (!aodTrack) AliFatal("Not a standard AOD");
2708 
2709  if (!aodTrack) continue;
2710 
2711  if (fAODTPCOnlyTracks)
2712  { // Match with TPC only tracks, default from May 2013, before filter bit 32
2713  //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
2714  if (!aodTrack->IsTPCConstrained()) continue ;
2715  }
2716  else if (fAODHybridTracks)
2717  { // Match with hybrid tracks
2718  //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
2719  if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
2720  } else
2721  { // Match with tracks on a mask
2722  //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask));
2723  if (!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks
2724  }
2725 
2726  if (aodTrack->Pt() < fCutMinTrackPt) continue;
2727 
2728  if ( TMath::Abs(aodTrack->Eta()) > 0.9 ) continue;
2729 
2730  // Save some time and memory in case of no DCal present
2731  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2732  {
2733  Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
2734  if ( phi <= 10 || phi >= 250 ) continue;
2735  }
2736 
2737  Double_t pos[3],mom[3];
2738  aodTrack->GetXYZ(pos);
2739  aodTrack->GetPxPyPz(mom);
2740  AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",
2741  itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
2742 
2743  trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
2744 
2745  mcLabel = TMath::Abs(aodTrack->GetLabel());
2746  }
2747 
2748  //Return if the input data is not "AOD" or "ESD"
2749  else
2750  {
2751  AliWarning("Wrong input data type! Should be \"AOD\" or \"ESD\" ");
2752  if (clusterArray)
2753  {
2754  clusterArray->Clear();
2755  delete clusterArray;
2756  }
2757  return;
2758  }
2759 
2760  if (!trackParam) continue;
2761 
2762  //
2763  // Check if track comes from a particular MC generator, do not include it
2764  // if it is not a selected one
2765  //
2766  if( mc && fMCGenerToAcceptForTrack && fNMCGenerToAccept > 0 )
2767  {
2768  mc->GetCocktailGenerator(mcLabel,genName);
2769 
2770  Bool_t generOK = kFALSE;
2771  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
2772  {
2773  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
2774  }
2775 
2776  if ( !generOK ) continue;
2777  }
2778 
2779  // Extrapolate the track to EMCal surface, see AliEMCALRecoUtilsBase
2780  AliExternalTrackParam emcalParam(*trackParam);
2781  Float_t eta, phi, pt;
2782  if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2783  {
2784  if (aodevent && trackParam) delete trackParam;
2785  if (fITSTrackSA && trackParam) delete trackParam;
2786  continue;
2787  }
2788 
2789  if ( TMath::Abs(eta) > 0.75 )
2790  {
2791  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2792  continue;
2793  }
2794 
2795  // Save some time and memory in case of no DCal present
2796  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2797  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2798  {
2799  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2800  continue;
2801  }
2802 
2803  //Find matched clusters
2804  Int_t index = -1;
2805  Float_t dEta = -999, dPhi = -999;
2806  if (!clusterArr)
2807  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
2808  else
2809  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr , dEta, dPhi);
2810 
2811 
2812  if (index>-1)
2813  {
2814  fMatchedTrackIndex ->AddAt(itr,matched);
2815  fMatchedClusterIndex ->AddAt(index,matched);
2816  fResidualEta ->AddAt(dEta,matched);
2817  fResidualPhi ->AddAt(dPhi,matched);
2818  matched++;
2819  }
2820 
2821  if (aodevent && trackParam) delete trackParam;
2822  if (fITSTrackSA && trackParam) delete trackParam;
2823  }//track loop
2824 
2825  if (clusterArray)
2826  {
2827  clusterArray->Clear();
2828  delete clusterArray;
2829  }
2830 
2831  AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
2832 
2833  fMatchedTrackIndex ->Set(matched);
2834  fMatchedClusterIndex ->Set(matched);
2835  fResidualPhi ->Set(matched);
2836  fResidualEta ->Set(matched);
2837 }
2838 
2850 //________________________________________________________________________________
2852  const AliVEvent *event,
2853  const AliEMCALGeometry *geom,
2854  Float_t &dEta, Float_t &dPhi)
2855 {
2856  Int_t index = -1;
2857 
2858  if ( TMath::Abs(track->Eta()) > 0.9 ) return index;
2859 
2860  // Save some time and memory in case of no DCal present
2861  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2862  {
2863  Double_t phiV = track->Phi()*TMath::RadToDeg();
2864  if ( phiV <= 10 || phiV >= 250 ) return index;
2865  }
2866 
2867  AliExternalTrackParam *trackParam = 0;
2868  if (!fITSTrackSA)
2869  trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam()); // If TPC
2870  else
2871  trackParam = new AliExternalTrackParam(*track);
2872 
2873  if (!trackParam) return index;
2874  AliExternalTrackParam emcalParam(*trackParam);
2875 
2876  Float_t eta, phi, pt;
2877  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2878  {
2879  if (fITSTrackSA) delete trackParam;
2880  return index;
2881  }
2882 
2883  if ( TMath::Abs(eta) > 0.75 )
2884  {
2885  if (fITSTrackSA) delete trackParam;
2886  return index;
2887  }
2888 
2889  // Save some time and memory in case of no DCal present
2890  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2891  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2892  {
2893  if (fITSTrackSA) delete trackParam;
2894  return index;
2895  }
2896 
2897  TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
2898 
2899  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2900  {
2901  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2902  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
2903  clusterArr->AddAt(cluster,icl);
2904  }
2905 
2906  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
2907  clusterArr->Clear();
2908  delete clusterArr;
2909  if (fITSTrackSA) delete trackParam;
2910 
2911  return index;
2912 }
2913 
2924 //_______________________________________________________________________________________________
2925 Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
2926  AliExternalTrackParam *trkParam,
2927  const TObjArray * clusterArr,
2928  Float_t &dEta, Float_t &dPhi)
2929 {
2930  dEta=-999, dPhi=-999;
2931  Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
2932  Int_t index = -1;
2933  Float_t tmpEta=-999, tmpPhi=-999;
2934 
2935  Double_t exPos[3] = {0.,0.,0.};
2936  if (!emcalParam->GetXYZ(exPos)) return index;
2937 
2938  Float_t clsPos[3] = {0.,0.,0.};
2939  for (Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
2940  {
2941  AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
2942 
2943  if (!cluster || !cluster->IsEMCAL()) continue;
2944 
2945  cluster->GetPosition(clsPos);
2946 
2947  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));
2948  if (dR > fClusterWindow) continue;
2949 
2950  AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
2951 
2952  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
2953 
2954  if (fCutEtaPhiSum)
2955  {
2956  Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
2957  if (tmpR<dRMax)
2958  {
2959  dRMax=tmpR;
2960  dEtaMax=tmpEta;
2961  dPhiMax=tmpPhi;
2962  index=icl;
2963  }
2964  }
2965  else if (fCutEtaPhiSeparate)
2966  {
2967  if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
2968  {
2969  dEtaMax = tmpEta;
2970  dPhiMax = tmpPhi;
2971  index=icl;
2972  }
2973  }
2974  else
2975  {
2976  AliError("Please specify your cut criteria");
2977  AliError("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()");
2978  AliError("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()");
2979  return index;
2980  }
2981  }
2982 
2983  dEta=dEtaMax;
2984  dPhi=dPhiMax;
2985 
2986  return index;
2987 }
2988 
3001 //---------------------------------------------------------------------------------
3002 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
3003  const AliVCluster *cluster,
3004  Float_t &tmpEta,
3005  Float_t &tmpPhi)
3006 {
3007  return AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
3008 }
3009 
3019 //_______________________________________________________________________
3021  Float_t &dEta, Float_t &dPhi)
3022 {
3023  if (FindMatchedPosForCluster(clsIndex) >= 999)
3024  {
3025  AliDebug(2,"No matched tracks found!\n");
3026  dEta=999.;
3027  dPhi=999.;
3028  return;
3029  }
3030 
3031  dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
3032  dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
3033 }
3034 
3044 //______________________________________________________________________________________________
3046 {
3047  if (FindMatchedPosForTrack(trkIndex) >= 999)
3048  {
3049  AliDebug(2,"No matched cluster found!\n");
3050  dEta=999.;
3051  dPhi=999.;
3052  return;
3053  }
3054 
3055  dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
3056  dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
3057 }
3058 
3065 //__________________________________________________________
3067 {
3068  if (IsClusterMatched(clsIndex))
3069  return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
3070  else
3071  return -1;
3072 }
3073 
3080 //__________________________________________________________
3082 {
3083  if (IsTrackMatched(trkIndex))
3084  return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
3085  else
3086  return -1;
3087 }
3088 
3096 //______________________________________________________________
3098 {
3099  if (FindMatchedPosForCluster(clsIndex) < 999)
3100  return kTRUE;
3101  else
3102  return kFALSE;
3103 }
3104 
3112 //____________________________________________________________
3114 {
3115  if (FindMatchedPosForTrack(trkIndex) < 999)
3116  return kTRUE;
3117  else
3118  return kFALSE;
3119 }
3120 
3128 //______________________________________________________________________
3130 {
3131  Float_t tmpR = fCutR;
3132  UInt_t pos = 999;
3133 
3134  for (Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
3135  {
3136  if (fMatchedClusterIndex->At(i)==clsIndex)
3137  {
3138  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3139  if (r<tmpR)
3140  {
3141  pos=i;
3142  tmpR=r;
3143  AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3144  fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3145  }
3146  }
3147  }
3148  return pos;
3149 }
3150 
3158 //____________________________________________________________________
3160 {
3161  Float_t tmpR = fCutR;
3162  UInt_t pos = 999;
3163 
3164  for (Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
3165  {
3166  if (fMatchedTrackIndex->At(i)==trkIndex)
3167  {
3168  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3169  if (r<tmpR)
3170  {
3171  pos=i;
3172  tmpR=r;
3173  AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3174  fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3175  }
3176  }
3177  }
3178  return pos;
3179 }
3180 
3195 //__________________________________________________________________________
3197  const AliEMCALGeometry *geom,
3198  AliVCaloCells* cells, Int_t bc)
3199 {
3200  Bool_t isGood=kTRUE;
3201 
3202  if (!cluster || !cluster->IsEMCAL()) return kFALSE;
3203  if (ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
3204  if (!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
3205  if (IsExoticCluster(cluster, cells,bc)) return kFALSE;
3206 
3207  return isGood;
3208 }
3209 
3220 //__________________________________________________________
3222 {
3223  UInt_t status = esdTrack->GetStatus();
3224 
3225  Int_t nClustersITS = esdTrack->GetITSclusters(0);
3226  Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
3227 
3228  Float_t chi2PerClusterITS = -1;
3229  Float_t chi2PerClusterTPC = -1;
3230  if (nClustersITS!=0)
3231  chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
3232  if (nClustersTPC!=0)
3233  chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
3234 
3235  //
3236  // DCA cuts
3237  // Only to be used for primary
3238  //
3239  if ( fTrackCutsType == kGlobalCut )
3240  {
3241  Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01);
3242  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
3243 
3244  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3245 
3246  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3247  }
3248  else if( fTrackCutsType == kGlobalCut2011 )
3249  {
3250  Float_t maxDCAToVertexXYPtDep = 0.0105 + 0.0350/TMath::Power(esdTrack->Pt(),1.1);
3251  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2011()
3252 
3253  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3254 
3255  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3256  }
3257 
3258  Float_t b[2];
3259  Float_t bCov[3];
3260  esdTrack->GetImpactParameters(b,bCov);
3261  if (bCov[0]<=0 || bCov[2]<=0)
3262  {
3263  AliDebug(1, "Estimated b resolution lower or equal zero!");
3264  bCov[0]=0; bCov[2]=0;
3265  }
3266 
3267  Float_t dcaToVertexXY = b[0];
3268  Float_t dcaToVertexZ = b[1];
3269  Float_t dcaToVertex = -1;
3270 
3271  if (fCutDCAToVertex2D)
3272  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY / fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY +
3273  dcaToVertexZ*dcaToVertexZ / fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ );
3274  else
3275  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
3276 
3277  // cut the track?
3278 
3279  Bool_t cuts[kNCuts];
3280  for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
3281 
3282  // track quality cuts
3283  if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
3284  cuts[0]=kTRUE;
3285  if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
3286  cuts[1]=kTRUE;
3287  if (nClustersTPC<fCutMinNClusterTPC)
3288  cuts[2]=kTRUE;
3289  if (nClustersITS<fCutMinNClusterITS)
3290  cuts[3]=kTRUE;
3291  if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
3292  cuts[4]=kTRUE;
3293  if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
3294  cuts[5]=kTRUE;
3295  if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
3296  cuts[6]=kTRUE;
3297  if (fCutDCAToVertex2D && dcaToVertex > 1)
3298  cuts[7] = kTRUE;
3299  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
3300  cuts[8] = kTRUE;
3301  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
3302  cuts[9] = kTRUE;
3303 
3305  {
3306  //Require at least one SPD point + anything else in ITS
3307  if ( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
3308  cuts[10] = kTRUE;
3309  }
3310 
3311  // ITS
3313  {
3314  if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin))
3315  {
3316  // TPC tracks
3317  cuts[11] = kTRUE;
3318  }
3319  else
3320  {
3321  // ITS standalone tracks
3323  {
3324  if (status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
3325  }
3326  else if (fCutRequireITSpureSA)
3327  {
3328  if (!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
3329  }
3330  }
3331  }
3332 
3333  Bool_t cut=kFALSE;
3334  for (Int_t i=0; i<kNCuts; i++)
3335  if (cuts[i]) { cut = kTRUE ; }
3336 
3337  // cut the track
3338  if (cut)
3339  return kFALSE;
3340  else
3341  return kTRUE;
3342 }
3343 
3349 //_____________________________________
3351 {
3352  switch (fTrackCutsType)
3353  {
3354  case kTPCOnlyCut:
3355  {
3356  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()"));
3357  //TPC
3358  SetMinNClustersTPC(70);
3360  SetAcceptKinkDaughters(kFALSE);
3361  SetRequireTPCRefit(kFALSE);
3362 
3363  //ITS
3364  SetRequireITSRefit(kFALSE);
3365  SetMaxDCAToVertexZ(3.2);
3366  SetMaxDCAToVertexXY(2.4);
3367  SetDCAToVertex2D(kTRUE);
3368 
3369  break;
3370  }
3371 
3372  case kGlobalCut:
3373  {
3374  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE)"));
3375  //TPC
3376  SetMinNClustersTPC(70);
3378  SetAcceptKinkDaughters(kFALSE);
3379  SetRequireTPCRefit(kTRUE);
3380 
3381  //ITS
3382  SetRequireITSRefit(kTRUE);
3383  SetMaxDCAToVertexZ(2);
3385  SetDCAToVertex2D(kFALSE);
3386 
3387  break;
3388  }
3389 
3390  case kLooseCut:
3391  {
3392  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
3393  SetMinNClustersTPC(50);
3394  SetAcceptKinkDaughters(kTRUE);
3395 
3396  break;
3397  }
3398 
3399  case kITSStandAlone:
3400  {
3401  AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
3402  SetRequireITSRefit(kTRUE);
3403  SetRequireITSStandAlone(kTRUE);
3404  SetITSTrackSA(kTRUE);
3405  break;
3406  }
3407 
3408  case kGlobalCut2011:
3409  {
3410  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE)"));
3411  //TPC
3412  SetMinNClustersTPC(50);
3414  SetAcceptKinkDaughters(kFALSE);
3415  SetRequireTPCRefit(kTRUE);
3416 
3417  //ITS
3418  SetRequireITSRefit(kTRUE);
3419  SetMaxDCAToVertexZ(2);
3421  SetDCAToVertex2D(kFALSE);
3422 
3423  break;
3424  }
3425 
3426  case kLooseCutWithITSrefit:
3427  {
3428  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut plus ITSrefit"));
3429  SetMinNClustersTPC(50);
3430  SetAcceptKinkDaughters(kTRUE);
3431  SetRequireITSRefit(kTRUE);
3432 
3433  break;
3434  }
3435  }
3436 }
3437 
3443 //________________________________________________________________________
3445 {
3446  Int_t nTracks = event->GetNumberOfTracks();
3447  for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
3448  {
3449  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
3450  if (!track)
3451  {
3452  AliWarning(Form("Could not receive track %d", iTrack));
3453  continue;
3454  }
3455 
3456  Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
3457  track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
3458 
3459  /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
3460  AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
3461  if (esdtrack)
3462  {
3463  if (matchClusIndex != -1)
3464  esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
3465  else
3466  esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3467  }
3468  else
3469  {
3470  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
3471  if (matchClusIndex != -1)
3472  aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
3473  else
3474  aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3475  }
3476  }
3477  AliDebug(2,"Track matched to closest cluster");
3478 }
3479 
3486 //_________________________________________________________________________
3488 {
3489  for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
3490  {
3491  AliVCluster *cluster = event->GetCaloCluster(iClus);
3492  if (!cluster->IsEMCAL())
3493  continue;
3494 
3495  //
3496  // Remove old matches in cluster
3497  //
3498  if(cluster->GetNTracksMatched() > 0)
3499  {
3500  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
3501  {
3502  TArrayI arrayTrackMatched(0);
3503  ((AliESDCaloCluster*)cluster)->AddTracksMatched(arrayTrackMatched);
3504  }
3505  else
3506  {
3507  for(Int_t iTrack = 0; iTrack < cluster->GetNTracksMatched(); iTrack++)
3508  {
3509  ((AliAODCaloCluster*)cluster)->RemoveTrackMatched((TObject*)((AliAODCaloCluster*)cluster)->GetTrackMatched(iTrack));
3510  }
3511  }
3512  }
3513 
3514  //
3515  // Find new matches and put them in the cluster
3516  //
3517  Int_t nTracks = event->GetNumberOfTracks();
3518  TArrayI arrayTrackMatched(nTracks);
3519 
3520  // Get the closest track matched to the cluster
3521  Int_t nMatched = 0;
3522  Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
3523  if (matchTrackIndex != -1)
3524  {
3525  arrayTrackMatched[nMatched] = matchTrackIndex;
3526  nMatched++;
3527  }
3528 
3529  // Get all other tracks matched to the cluster
3530  for (Int_t iTrk=0; iTrk<nTracks; ++iTrk)
3531  {
3532  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
3533 
3534  if( !track ) continue;
3535 
3536  if ( iTrk == matchTrackIndex ) continue;
3537 
3538  if ( track->GetEMCALcluster() == iClus )
3539  {
3540  arrayTrackMatched[nMatched] = iTrk;
3541  ++nMatched;
3542  }
3543  }
3544 
3545  AliDebug(2,Form("cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]));
3546 
3547  arrayTrackMatched.Set(nMatched);
3548  AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
3549  if (esdcluster)
3550  esdcluster->AddTracksMatched(arrayTrackMatched);
3551  else if ( nMatched > 0 )
3552  {
3553  AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
3554  if (aodcluster)
3555  {
3556  aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
3557  //AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(event->GetTrack(arrayTrackMatched.At(0)));
3558  //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());
3559  //printf("With specs: pt %.4f, eta %.4f, phi %.4f\n",aodtrack->Pt(),aodtrack->Eta(), aodtrack->Phi());
3560  }
3561  }
3562 
3563  Float_t eta= -999, phi = -999;
3564  if (matchTrackIndex != -1)
3565  GetMatchedResiduals(iClus, eta, phi);
3566 
3567  cluster->SetTrackDistance(phi, eta);
3568  }
3569 
3570  AliDebug(2,"Cluster matched to tracks");
3571 }
3572 
3575  else {
3576  fEMCALRecalibrationFactors = new TObjArray(map->GetEntries());
3577  fEMCALRecalibrationFactors->SetOwner(true);
3578  }
3579  if(!fEMCALRecalibrationFactors->IsOwner()){
3580  // Must claim ownership since the new objects are owend by this instance
3581  fEMCALRecalibrationFactors->SetOwner(kTRUE);
3582  }
3583  for(int i = 0; i < map->GetEntries(); i++){
3584  TH2F *hist = dynamic_cast<TH2F *>(map->At(i));
3585  if(!hist) continue;
3586  this->SetEMCALChannelRecalibrationFactors(i, hist);
3587  }
3588 }
3589 
3593  fEMCALRecalibrationFactors->SetOwner(true);
3594  }
3595  if(fEMCALRecalibrationFactors->GetEntries() <= iSM) fEMCALRecalibrationFactors->Expand(iSM+1);
3596  if(fEMCALRecalibrationFactors->At(iSM)) fEMCALRecalibrationFactors->RemoveAt(iSM);
3597  TH2F *clone = new TH2F(*h);
3598  clone->SetDirectory(NULL);
3599  fEMCALRecalibrationFactors->AddAt(clone,iSM);
3600 }
3601 
3604  else {
3605  fEMCALBadChannelMap = new TObjArray(map->GetEntries());
3606  fEMCALBadChannelMap->SetOwner(true);
3607  }
3608  if(!fEMCALBadChannelMap->IsOwner()) {
3609  // Must claim ownership since the new objects are owend by this instance
3610  fEMCALBadChannelMap->SetOwner(true);
3611  }
3612  for(int i = 0; i < map->GetEntries(); i++){
3613  TH2I *hist = dynamic_cast<TH2I *>(map->At(i));
3614  if(!hist) continue;
3615  this->SetEMCALChannelStatusMap(i, hist);
3616  }
3617 }
3618 
3620  if(!fEMCALBadChannelMap){
3621  fEMCALBadChannelMap = new TObjArray(iSM);
3622  fEMCALBadChannelMap->SetOwner(true);
3623  }
3624  if(fEMCALBadChannelMap->GetEntries() <= iSM) fEMCALBadChannelMap->Expand(iSM+1);
3625  if(fEMCALBadChannelMap->At(iSM)) fEMCALBadChannelMap->RemoveAt(iSM);
3626  TH2I *clone = new TH2I(*h);
3627  clone->SetDirectory(NULL);
3628  fEMCALBadChannelMap->AddAt(clone,iSM);
3629 }
3630 
3633  else {
3634  fEMCALTimeRecalibrationFactors = new TObjArray(map->GetEntries());
3635  fEMCALTimeRecalibrationFactors->SetOwner(true);
3636  }
3637  if(!fEMCALTimeRecalibrationFactors->IsOwner()) {
3638  // Must claim ownership since the new objects are owend by this instance
3639  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
3640  }
3641  for(int i = 0; i < map->GetEntries(); i++){
3642  TH1F *hist = dynamic_cast<TH1F *>(map->At(i));
3643  if(!hist) continue;
3645  }
3646 }
3647 
3651  fEMCALTimeRecalibrationFactors->SetOwner(true);
3652  }
3653  if(fEMCALTimeRecalibrationFactors->GetEntries() <= bc) fEMCALTimeRecalibrationFactors->Expand(bc+1);
3655  TH1F *clone = new TH1F(*h);
3656  clone->SetDirectory(NULL);
3657  fEMCALTimeRecalibrationFactors->AddAt(clone,bc);
3658 }
3659 
3662  else {
3663  fEMCALL1PhaseInTimeRecalibration = new TObjArray(map->GetEntries());
3664  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3665  }
3666  if(!fEMCALL1PhaseInTimeRecalibration->IsOwner()) {
3667  // Must claim ownership since the new objects are owend by this instance
3668  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3669  }
3670  for(int i = 0; i < map->GetEntries(); i++) {
3671  TH1C *hist = dynamic_cast<TH1C *>(map->At(i));
3672  if(!hist) continue;
3674  }
3675 }
3676 
3680  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3681  }
3683  TH1C *clone = new TH1C(*h);
3684  clone->SetDirectory(NULL);
3685  fEMCALL1PhaseInTimeRecalibration->AddAt(clone,0);
3686 }
3687 
3691 //___________________________________________________
3693 {
3694  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3695  printf("AliEMCALRecoUtils Settings: \n");
3696  printf("\tMisalignment shifts\n");
3697  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,
3699  fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
3700  printf("\tNon linearity function %d, parameters:\n", fNonLinearityFunction);
3701  if (fNonLinearityFunction != 3) // print only if not kNoCorrection
3702  for (Int_t i=0; i<10; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
3703 
3704  printf("\tPosition Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
3705 
3706  printf("\tMatching criteria: ");
3707  if (fCutEtaPhiSum) {
3708  printf("\t\tsqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
3709  } else if (fCutEtaPhiSeparate) {
3710  printf("\t\tdEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
3711  } else {
3712  printf("\t\tError\n");
3713  printf("\t\tplease specify your cut criteria\n");
3714  printf("\t\tTo cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
3715  printf("\t\tTo cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
3716  }
3717 
3718  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);
3719  printf("\tCluster selection window: dR < %2.0f\n",fClusterWindow);
3720 
3721  printf("\tTrack cuts: \n");
3722  printf("\t\tMinimum track pT: %1.2f\n",fCutMinTrackPt);
3723  printf("\t\tAOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask);
3724  printf("\t\tTPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
3725  printf("\t\tAcceptKinks = %d\n",fCutAcceptKinkDaughters);
3726  printf("\t\tMinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
3727  printf("\t\tMaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
3728  printf("\t\tDCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
3729  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3730 }
void RecalculateCellLabelsRemoveAddedGenerator(Int_t absID, AliVCluster *clus, AliMCEvent *mc, Float_t &amp, TArrayI &labeArr, TArrayF &eDepArr) const
Bool_t IsL1PhaseInTimeRecalibrationOn() const
Double_t fMass
Mass hypothesis of the track.
Bool_t fMCGenerToAcceptForTrack
Activate the removal of tracks entering the track matching that come from a particular generator...
double Double_t
Definition: External.C:58
Bool_t CheckCellFiducialRegion(const AliEMCALGeometry *geom, const AliVCluster *cluster, AliVCaloCells *cells)
Bool_t fAODTPCOnlyTracks
Match with TPC only tracks.
Bool_t fCutRequireTPCRefit
Require TPC refit.
Definition: External.C:236
void SetRequireITSStandAlone(Bool_t b=kFALSE)
Bool_t IsClusterMatched(Int_t clsIndex) const
Float_t GetEMCALChannelTimeRecalibrationFactor(Int_t bc, Int_t absID, Bool_t isLGon=kFALSE) const
TArrayI * fMatchedClusterIndex
Array that stores indexes of matched clusters.
TObjArray * fEMCALTimeRecalibrationFactors
Array of histograms with map of time recalibration factors, EMCAL.
void SetMaxDCAToVertexXY(Float_t dist=1e10)
TRandom3 fRandom
Random generator for cluster energy smearing.
Bool_t fTimeRecalibration
Switch on or off the time recalibration.
Bool_t fCutEtaPhiSum
Place cut on sqrt(dEta^2+dPhi^2)
Int_t FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam, AliExternalTrackParam *trkParam, const TObjArray *clusterArr, Float_t &dEta, Float_t &dPhi)
void InitEMCALL1PhaseInTimeRecalibration()
Int_t fCutMinNClusterTPC
Min number of tpc clusters.
Float_t fCutMaxDCAToVertexZ
Track-to-vertex cut in max absolute distance in z-plane.
void SetRequireITSRefit(Bool_t b=kFALSE)
energy
Definition: HFPtSpectrum.C:44
void SetEMCALL1PhaseInTimeRecalibrationForSM(Int_t iSM, Int_t c=0)
void SetDCAToVertex2D(Bool_t b=kFALSE)
void SetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow, Double_t c=1)
Int_t fNonLinearThreshold
Non linearity threshold value for kBeamTest non linearity function.
Bool_t fITSTrackSA
If track matching is to be done with ITS tracks standing alone.
void SetITSTrackSA(Bool_t isITS)
Bool_t IsExoticCluster(const AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=0)
Float_t fMisalRotShift[15]
Cluster position rotation shift parameters.
void SetMinNClustersTPC(Int_t min=-1)
Bool_t fRecalDistToBadChannels
Calculate distance from highest energy tower of cluster to closes bad channel.
Bool_t fRecalibration
Switch on or off the recalibration.
TArrayI * fMatchedTrackIndex
Array that stores indexes of matched tracks.
void SetMaxDCAToVertexZ(Float_t dist=1e10)
Bool_t IsBadChannelsRemovalSwitchedOn() const
Bool_t IsExoticCell(Int_t absId, AliVCaloCells *cells, Int_t bc=-1)
void RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
Float_t fExoticCellDiffTime
If time of candidate to exotic and close cell is too different (in ns), it must be noisy...
TArrayF * fResidualPhi
Array that stores the residual phi.
void RecalibrateCells(AliVCaloCells *cells, Int_t bc)
Double_t fCutMinTrackPt
Cut on track pT.
Bool_t IsRecalibrationOn() const
Some utilities for cluster and cell treatment.
Float_t fCutMaxChi2PerClusterITS
Max its fit chi2 per its cluster.
void RecalculateClusterShowerShapeParametersWithCellCuts(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster, Float_t cellEcut, Float_t cellTimeCut, Int_t bc, Float_t &enAfterCuts)
Bool_t fRemoveBadChannels
Check the channel status provided and remove clusters with bad channels.
Bool_t IsAccepted(AliESDtrack *track)
Float_t fMisalTransShift[15]
Cluster position translation shift parameters.
void SetAcceptKinkDaughters(Bool_t b=kTRUE)
void RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
Double_t fClusterWindow
Select clusters in the window to be matched.
Bool_t fCutRequireITSRefit
Require ITS refit.
Float_t fNonLinearityParams[10]
Parameters for the non linearity function.
Float_t GetCellWeight(Float_t eCell, Float_t eCluster) const
Bool_t fUseL1PhaseInTimeRecalibration
Switch on or off the L1 phase in time recalibration.
Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow) const
void SetClusterMatchedToTrack(const AliVEvent *event)
Bool_t fCutDCAToVertex2D
If true a 2D DCA cut is made.
Int_t GetMatchedTrackIndex(Int_t clsIndex)
void RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
int Int_t
Definition: External.C:63
void Print(const Option_t *) const
UInt_t FindMatchedPosForCluster(Int_t clsIndex) const
Bool_t AcceptCalibrateCell(Int_t absId, Int_t bc, Float_t &amp, Double_t &time, AliVCaloCells *cells)
unsigned int UInt_t
Definition: External.C:33
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
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)