AliPhysics  b76e98e (b76e98e)
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 parametrization 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 kBeamTestCorrected, 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 
916  {
917  // New parametrization of kBeamTestCorrected,
918  // fitting new points for E>100 GeV.
919  // I should have same performance as v3 in the low energies
920  // See EMCal meeting 21/09/2018 slides
921  // https://indico.cern.ch/event/759154/contributions/3148448/attachments/1721042/2778585/nonLinearityUpdate.pdf
922  // and jira ticket EMCAL-190
923  // Not very smart copy pasting the same function for each new parametrization, need to think how to do it better.
924 
925 // fNonLinearityParams[0] = 0.9892;
926 // fNonLinearityParams[1] = 0.1976;
927 // fNonLinearityParams[2] = 0.865;
928 // fNonLinearityParams[3] = 0.06775;
929 // fNonLinearityParams[4] = 156.6;
930 // fNonLinearityParams[5] = 47.18;
931 // fNonLinearityParams[6] = 0.97;
932 
933  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
934 
935  break;
936  }
937 
938  case kSDMv5:
939  {
940  //Based on fit to the MC/data using kNoCorrection on the data - utilizes symmetric decay method and kPi0MCv5(MC) - 28 Oct 2013
941  //fNonLinearityParams[0] = 1.0;
942  //fNonLinearityParams[1] = 6.64778e-02;
943  //fNonLinearityParams[2] = 1.570;
944  //fNonLinearityParams[3] = 9.67998e-02;
945  //fNonLinearityParams[4] = 2.19381e+02;
946  //fNonLinearityParams[5] = 6.31604e+01;
947  //fNonLinearityParams[6] = 1.01286;
948  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));
949 
950  break;
951  }
952 
953  case kPi0MCv5:
954  {
955  //Based on comparing MC truth information to the reconstructed energy of clusters.
956  //fNonLinearityParams[0] = 1.0;
957  //fNonLinearityParams[1] = 6.64778e-02;
958  //fNonLinearityParams[2] = 1.570;
959  //fNonLinearityParams[3] = 9.67998e-02;
960  //fNonLinearityParams[4] = 2.19381e+02;
961  //fNonLinearityParams[5] = 6.31604e+01;
962  //fNonLinearityParams[6] = 1.01286;
963  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
964 
965  break;
966  }
967 
968  case kSDMv6:
969  {
970  //Based on fit to the MC/data using kNoCorrection on the data
971  // - utilizes symmetric decay method and kPi0MCv6(MC) - 09 Dec 2014
972  // - parameters constrained by the test beam data as well
973  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
974  //fNonLinearityParams[0] = 1.0;
975  //fNonLinearityParams[1] = 0.237767;
976  //fNonLinearityParams[2] = 0.651203;
977  //fNonLinearityParams[3] = 0.183741;
978  //fNonLinearityParams[4] = 155.427;
979  //fNonLinearityParams[5] = 17.0335;
980  //fNonLinearityParams[6] = 0.987054;
981  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
982 
983  break;
984  }
985 
986  case kPi0MCv6:
987  {
988  //Based on comparing MC truth information to the reconstructed energy of clusters.
989  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
990  //fNonLinearityParams[0] = 1.0;
991  //fNonLinearityParams[1] = 0.0797873;
992  //fNonLinearityParams[2] = 1.68322;
993  //fNonLinearityParams[3] = 0.0806098;
994  //fNonLinearityParams[4] = 244.586;
995  //fNonLinearityParams[5] = 116.938;
996  //fNonLinearityParams[6] = 1.00437;
997  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
998 
999  break;
1000  }
1001 
1002  case kPCMv1:
1003  {
1004  //based on symmetric decays of pi0 meson
1005  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1006  // parameters vary from MC to MC
1007  //fNonLinearityParams[0] = 0.984876;
1008  //fNonLinearityParams[1] = -9.999609;
1009  //fNonLinearityParams[2] = -4.999891;
1010  //fNonLinearityParams[3] = 0.;
1011  //fNonLinearityParams[4] = 0.;
1012  //fNonLinearityParams[5] = 0.;
1013  //fNonLinearityParams[6] = 0.;
1014  energy /= TMath::Power(fNonLinearityParams[0] + exp(fNonLinearityParams[1] + fNonLinearityParams[2]*energy),2);//result coming from calo-conv method
1015 
1016  break;
1017  }
1018 
1019  case kPCMplusBTCv1:
1020  {
1021  //convolution of TestBeamCorrectedv3 with PCM method
1022  //Based on comparing MC truth information to the reconstructed energy of clusters.
1023  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1024  // parameters vary from MC to MC
1025  //fNonLinearityParams[0] = 0.976941;
1026  //fNonLinearityParams[1] = 0.162310;
1027  //fNonLinearityParams[2] = 1.08689;
1028  //fNonLinearityParams[3] = 0.0819592;
1029  //fNonLinearityParams[4] = 152.338;
1030  //fNonLinearityParams[5] = 30.9594;
1031  //fNonLinearityParams[6] = 0.9615;
1032  //fNonLinearityParams[7] = 0.984876;
1033  //fNonLinearityParams[8] = -9.999609;
1034  //fNonLinearityParams[9] = -4.999891;
1035  energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
1036  energy /= TMath::Power(fNonLinearityParams[7] + exp(fNonLinearityParams[8] + fNonLinearityParams[9]*energy),2);//result coming from calo-conv method
1037 
1038  break;
1039  }
1040 
1041  case kPCMsysv1:
1042  {
1043  // Systematic variation of kPCMv1
1044  //Based on comparing MC truth information to the reconstructed energy of clusters.
1045  // described in the note: https://aliceinfo.cern.ch/Notes/node/211 - Sec 3.1.2 (Test Beam Constrained SDM).
1046  // parameters vary from MC to MC
1047  //fNonLinearityParams[0] = 0.0;
1048  //fNonLinearityParams[1] = 1.0;
1049  //fNonLinearityParams[2] = 1.0;
1050  //fNonLinearityParams[3] = 0.0;
1051  //fNonLinearityParams[4] = 1.0;
1052  //fNonLinearityParams[5] = 0.0;
1053  //fNonLinearityParams[6] = 0.0;
1054  energy /= TMath::Power( (fNonLinearityParams[0] + fNonLinearityParams[1] * TMath::Power(energy,fNonLinearityParams[2]) ) /
1055  (fNonLinearityParams[3] + fNonLinearityParams[4] * TMath::Power(energy,fNonLinearityParams[5]) ) + fNonLinearityParams[6], 2);//result coming from calo-conv method
1056 
1057  break;
1058  }
1059 
1060 
1061  case kNoCorrection:
1062  AliDebug(2,"No correction on the energy\n");
1063  break;
1064 
1065  }
1066 
1067  return energy;
1068 }
1069 
1074 //__________________________________________________
1076 {
1077  if (fNonLinearityFunction == kPi0MC) {
1078  fNonLinearityParams[0] = 1.014;
1079  fNonLinearityParams[1] = -0.03329;
1080  fNonLinearityParams[2] = -0.3853;
1081  fNonLinearityParams[3] = 0.5423;
1082  fNonLinearityParams[4] = -0.4335;
1083  }
1084 
1086  fNonLinearityParams[0] = 3.11111e-02;
1087  fNonLinearityParams[1] =-5.71666e-02;
1088  fNonLinearityParams[2] = 5.67995e-01;
1089  }
1090 
1092  fNonLinearityParams[0] = 9.81039e-01;
1093  fNonLinearityParams[1] = 1.13508e-01;
1094  fNonLinearityParams[2] = 1.00173e+00;
1095  fNonLinearityParams[3] = 9.67998e-02;
1096  fNonLinearityParams[4] = 2.19381e+02;
1097  fNonLinearityParams[5] = 6.31604e+01;
1098  fNonLinearityParams[6] = 1;
1099  }
1100 
1102  fNonLinearityParams[0] = 1.04;
1103  fNonLinearityParams[1] = -0.1445;
1104  fNonLinearityParams[2] = 1.046;
1105  }
1106 
1108  fNonLinearityParams[0] = 0.139393;
1109  fNonLinearityParams[1] = 0.0566186;
1110  fNonLinearityParams[2] = 0.982133;
1111  }
1112 
1114  if (fNonLinearThreshold == 30) {
1115  fNonLinearityParams[0] = 1.007;
1116  fNonLinearityParams[1] = 0.894;
1117  fNonLinearityParams[2] = 0.246;
1118  }
1119  if (fNonLinearThreshold == 45) {
1120  fNonLinearityParams[0] = 1.003;
1121  fNonLinearityParams[1] = 0.719;
1122  fNonLinearityParams[2] = 0.334;
1123  }
1124  if (fNonLinearThreshold == 75) {
1125  fNonLinearityParams[0] = 1.002;
1126  fNonLinearityParams[1] = 0.797;
1127  fNonLinearityParams[2] = 0.358;
1128  }
1129  }
1130 
1132  fNonLinearityParams[0] = 0.99078;
1133  fNonLinearityParams[1] = 0.161499;
1134  fNonLinearityParams[2] = 0.655166;
1135  fNonLinearityParams[3] = 0.134101;
1136  fNonLinearityParams[4] = 163.282;
1137  fNonLinearityParams[5] = 23.6904;
1138  fNonLinearityParams[6] = 0.978;
1139  }
1140 
1142  // Parameters until November 2015, use now kBeamTestCorrectedv3
1143  fNonLinearityParams[0] = 0.983504;
1144  fNonLinearityParams[1] = 0.210106;
1145  fNonLinearityParams[2] = 0.897274;
1146  fNonLinearityParams[3] = 0.0829064;
1147  fNonLinearityParams[4] = 152.299;
1148  fNonLinearityParams[5] = 31.5028;
1149  fNonLinearityParams[6] = 0.968;
1150  }
1151 
1153 
1154  // New parametrization of kBeamTestCorrected
1155  // excluding point at 0.5 GeV from Beam Test Data
1156  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1157 
1158  fNonLinearityParams[0] = 0.976941;
1159  fNonLinearityParams[1] = 0.162310;
1160  fNonLinearityParams[2] = 1.08689;
1161  fNonLinearityParams[3] = 0.0819592;
1162  fNonLinearityParams[4] = 152.338;
1163  fNonLinearityParams[5] = 30.9594;
1164  fNonLinearityParams[6] = 0.9615;
1165  }
1166 
1168 
1169  // New parametrization of kBeamTestCorrected,
1170  // fitting new points for E>100 GeV.
1171  // I should have same performance as v3 in the low energies
1172  // See EMCal meeting 21/09/2018 slides
1173  // https://indico.cern.ch/event/759154/contributions/3148448/attachments/1721042/2778585/nonLinearityUpdate.pdf
1174  // and jira ticket EMCAL-190
1175 
1176  fNonLinearityParams[0] = 0.9892;
1177  fNonLinearityParams[1] = 0.1976;
1178  fNonLinearityParams[2] = 0.865;
1179  fNonLinearityParams[3] = 0.06775;
1180  fNonLinearityParams[4] = 156.6;
1181  fNonLinearityParams[5] = 47.18;
1182  fNonLinearityParams[6] = 0.97;
1183  }
1184 
1185  if (fNonLinearityFunction == kSDMv5) {
1186  fNonLinearityParams[0] = 1.0;
1187  fNonLinearityParams[1] = 6.64778e-02;
1188  fNonLinearityParams[2] = 1.570;
1189  fNonLinearityParams[3] = 9.67998e-02;
1190  fNonLinearityParams[4] = 2.19381e+02;
1191  fNonLinearityParams[5] = 6.31604e+01;
1192  fNonLinearityParams[6] = 1.01286;
1193  }
1194 
1196  fNonLinearityParams[0] = 1.0;
1197  fNonLinearityParams[1] = 6.64778e-02;
1198  fNonLinearityParams[2] = 1.570;
1199  fNonLinearityParams[3] = 9.67998e-02;
1200  fNonLinearityParams[4] = 2.19381e+02;
1201  fNonLinearityParams[5] = 6.31604e+01;
1202  fNonLinearityParams[6] = 1.01286;
1203  }
1204 
1205  if (fNonLinearityFunction == kSDMv6) {
1206  fNonLinearityParams[0] = 1.0;
1207  fNonLinearityParams[1] = 0.237767;
1208  fNonLinearityParams[2] = 0.651203;
1209  fNonLinearityParams[3] = 0.183741;
1210  fNonLinearityParams[4] = 155.427;
1211  fNonLinearityParams[5] = 17.0335;
1212  fNonLinearityParams[6] = 0.987054;
1213  }
1214 
1216  fNonLinearityParams[0] = 1.0;
1217  fNonLinearityParams[1] = 0.0797873;
1218  fNonLinearityParams[2] = 1.68322;
1219  fNonLinearityParams[3] = 0.0806098;
1220  fNonLinearityParams[4] = 244.586;
1221  fNonLinearityParams[5] = 116.938;
1222  fNonLinearityParams[6] = 1.00437;
1223  }
1224 
1225 if (fNonLinearityFunction == kPCMv1) {
1226  //parameters change from MC production to MC production, they need to set for each period
1227  fNonLinearityParams[0] = 0.984876;
1228  fNonLinearityParams[1] = -9.999609;
1229  fNonLinearityParams[2] = -4.999891;
1230  fNonLinearityParams[3] = 0.;
1231  fNonLinearityParams[4] = 0.;
1232  fNonLinearityParams[5] = 0.;
1233  fNonLinearityParams[6] = 0.;
1234  }
1235 
1237  // test beam corrected values convoluted with symmetric meson decays values
1238  // for test beam:
1239  // https://indico.cern.ch/event/438805/contribution/1/attachments/1145354/1641875/emcalPi027August2015.pdf
1240  // for PCM method:
1241  // https://aliceinfo.cern.ch/Notes/node/211
1242  fNonLinearityParams[0] = 0.976941;
1243  fNonLinearityParams[1] = 0.162310;
1244  fNonLinearityParams[2] = 1.08689;
1245  fNonLinearityParams[3] = 0.0819592;
1246  fNonLinearityParams[4] = 152.338;
1247  fNonLinearityParams[5] = 30.9594;
1248  fNonLinearityParams[6] = 0.9615;
1249  fNonLinearityParams[7] = 0.984876;
1250  fNonLinearityParams[8] = -9.999609;
1251  fNonLinearityParams[9] = -4.999891;
1252  }
1253 
1255  //systematics for kPCMv1
1256  // for PCM method:
1257  // https://aliceinfo.cern.ch/Notes/node/211
1258  fNonLinearityParams[0] = 0.0;
1259  fNonLinearityParams[1] = 1.0;
1260  fNonLinearityParams[2] = 1.0;
1261  fNonLinearityParams[3] = 0.0;
1262  fNonLinearityParams[4] = 1.0;
1263  fNonLinearityParams[5] = 0.0;
1264  fNonLinearityParams[6] = 0.0;
1265  }
1266 }
1267 
1275 //____________________________________________________________________
1277 {
1278  fBadStatusSelection[0] = all; // declare all as bad if true, never mind the other settings
1279  fBadStatusSelection[1] = dead;
1280  fBadStatusSelection[2] = hot;
1281  fBadStatusSelection[3] = warm;
1282 }
1283 
1294 //____________________________________________________________________
1296 {
1297  if(fEMCALBadChannelMap)
1298  status = (Int_t) ((TH2I*)fEMCALBadChannelMap->At(iSM))->GetBinContent(iCol,iRow);
1299  else
1300  status = 0; // Channel is ok by default
1301 
1302  if ( status == AliCaloCalibPedestal::kAlive )
1303  {
1304  return kFALSE; // Good channel
1305  }
1306  else
1307  {
1308  if ( fBadStatusSelection[0] == kTRUE )
1309  {
1310  return kTRUE; // consider bad hot, dead and warm
1311  }
1312  else
1313  {
1314  if ( fBadStatusSelection[AliCaloCalibPedestal::kDead] == kTRUE &&
1315  status == AliCaloCalibPedestal::kDead )
1316  return kTRUE; // consider bad dead
1317  else if ( fBadStatusSelection[AliCaloCalibPedestal::kHot] == kTRUE &&
1318  status == AliCaloCalibPedestal::kHot )
1319  return kTRUE; // consider bad hot
1320  else if ( fBadStatusSelection[AliCaloCalibPedestal::kWarning] == kTRUE &&
1321  status == AliCaloCalibPedestal::kWarning )
1322  return kTRUE; // consider bad warm
1323  }
1324  }
1325 
1326  AliWarning(Form("Careful, bad channel selection not properly done: ism %d, icol %d, irow %d, status %d,\n"
1327  " fBadAll %d, fBadHot %d, fBadWarm %d, fBadDead %d",
1328  iSM, iCol, iRow, status,
1331 
1332  return kFALSE; // if everything fails, accept it.
1333 }
1334 
1348 //____________________________________________________________________
1349 void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
1350  AliVCaloCells* cells,
1351  const AliVCluster* clu,
1352  Int_t & absId,
1353  Int_t & iSupMod,
1354  Int_t & ieta,
1355  Int_t & iphi,
1356  Bool_t & shared)
1357 {
1358  Double_t eMax = -1.;
1359  Double_t eCell = -1.;
1360  Float_t fraction = 1.;
1361  Float_t recalFactor = 1.;
1362  Int_t cellAbsId = -1 ;
1363 
1364  Int_t iTower = -1;
1365  Int_t iIphi = -1;
1366  Int_t iIeta = -1;
1367  Int_t iSupMod0= -1;
1368 
1369  if (!clu)
1370  {
1371  AliInfo("Cluster pointer null!");
1372  absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
1373  return;
1374  }
1375 
1376  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1377  {
1378  cellAbsId = clu->GetCellAbsId(iDig);
1379  fraction = clu->GetCellAmplitudeFraction(iDig);
1380  //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
1381  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
1382 
1383  geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1384  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1385 
1386  if (iDig==0)
1387  {
1388  iSupMod0=iSupMod;
1389  } else if (iSupMod0!=iSupMod)
1390  {
1391  shared = kTRUE;
1392  //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
1393  }
1394 
1396  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1397 
1398  eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
1399  //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
1400  if (eCell > eMax)
1401  {
1402  eMax = eCell;
1403  absId = cellAbsId;
1404  //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
1405  }
1406  }// cell loop
1407 
1408  //Get from the absid the supermodule, tower and eta/phi numbers
1409  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1410  //Gives SuperModule and Tower numbers
1411  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1412  //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
1413  //printf("Max end---\n");
1414 }
1415 
1424 //_________________________________________________________
1426 {
1427  if (eCell > 0 && eCluster > 0)
1428  {
1429  if ( fW0 > 0 ) return TMath::Max( 0., fW0 + TMath::Log( eCell / eCluster ) ) ;
1430  else return TMath::Log( eCluster / eCell ) ;
1431  }
1432  else
1433  return 0. ;
1434 }
1435 
1439 //______________________________________
1441 {
1442  fParticleType = kPhoton;
1443  fPosAlgo = kUnchanged;
1444  fW0 = 4.5;
1445 
1447  fNonLinearThreshold = 30;
1448 
1449  fExoticCellFraction = 0.97;
1450  fExoticCellDiffTime = 1e6;
1452 
1453  fAODFilterMask = 128;
1454  fAODHybridTracks = kFALSE;
1455  fAODTPCOnlyTracks = kTRUE;
1456 
1457  fCutEtaPhiSum = kTRUE;
1458  fCutEtaPhiSeparate = kFALSE;
1459 
1460  fCutR = 0.05;
1461  fCutEta = 0.025;
1462  fCutPhi = 0.05;
1463 
1464  fClusterWindow = 100;
1465  fMass = 0.139;
1466 
1467  fStepSurface = 20.;
1468  fStepCluster = 5.;
1470 
1471  fCutMinTrackPt = 0;
1472  fCutMinNClusterTPC = -1;
1473  fCutMinNClusterITS = -1;
1474 
1475  fCutMaxChi2PerClusterTPC = 1e10;
1476  fCutMaxChi2PerClusterITS = 1e10;
1477 
1478  fCutRequireTPCRefit = kFALSE;
1479  fCutRequireITSRefit = kFALSE;
1480  fCutAcceptKinkDaughters = kFALSE;
1481 
1482  fCutMaxDCAToVertexXY = 1e10;
1483  fCutMaxDCAToVertexZ = 1e10;
1484  fCutDCAToVertex2D = kFALSE;
1485 
1486  fCutRequireITSStandAlone = kFALSE; //MARCEL
1487  fCutRequireITSpureSA = kFALSE; //Marcel
1488 
1489  // Misalignment matrices
1490  for (Int_t i = 0; i < 15 ; i++)
1491  {
1492  fMisalTransShift[i] = 0.;
1493  fMisalRotShift[i] = 0.;
1494  }
1495 
1496  // Non linearity
1497  for (Int_t i = 0; i < 10 ; i++) fNonLinearityParams[i] = 0.;
1498 
1499  // For kBeamTestCorrectedv2 case, but default is no correction
1500  fNonLinearityParams[0] = 0.983504;
1501  fNonLinearityParams[1] = 0.210106;
1502  fNonLinearityParams[2] = 0.897274;
1503  fNonLinearityParams[3] = 0.0829064;
1504  fNonLinearityParams[4] = 152.299;
1505  fNonLinearityParams[5] = 31.5028;
1506  fNonLinearityParams[6] = 0.968;
1507 
1508  // Cluster energy smearing
1509  fSmearClusterEnergy = kFALSE;
1510  fSmearClusterParam[0] = 0.07; // * sqrt E term
1511  fSmearClusterParam[1] = 0.00; // * E term
1512  fSmearClusterParam[2] = 0.00; // constant
1513 }
1514 
1518 //_____________________________________________________
1520 {
1521  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1522 
1523  // In order to avoid rewriting the same histograms
1524  Bool_t oldStatus = TH1::AddDirectoryStatus();
1525  TH1::AddDirectory(kFALSE);
1526 
1528  for (int i = 0; i < 22; i++)
1529  fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1530  Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
1531  //Init the histograms with 1
1532  for (Int_t sm = 0; sm < 22; sm++)
1533  {
1534  for (Int_t i = 0; i < 48; i++)
1535  {
1536  for (Int_t j = 0; j < 24; j++)
1537  {
1539  }
1540  }
1541  }
1542 
1543  fEMCALRecalibrationFactors->SetOwner(kTRUE);
1544  fEMCALRecalibrationFactors->Compress();
1545 
1546  // In order to avoid rewriting the same histograms
1547  TH1::AddDirectory(oldStatus);
1548 }
1549 
1553 //_________________________________________________________
1555 {
1556  AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1557 
1558  // In order to avoid rewriting the same histograms
1559  Bool_t oldStatus = TH1::AddDirectoryStatus();
1560  TH1::AddDirectory(kFALSE);
1561 
1564 
1565  for (int i = 0; i < 4; i++)
1566  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1567  Form("hAllTimeAvBC%d",i),
1568  48*24*22,0.,48*24*22) );
1569  // Init the histograms with 0
1570  for (Int_t iBC = 0; iBC < 4; iBC++)
1571  {
1572  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1573  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kFALSE);
1574  }
1575 
1576  if(fLowGain) {
1577  for (int iBC = 0; iBC < 4; iBC++) {
1578  fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvLGBC%d",iBC),
1579  Form("hAllTimeAvLGBC%d",iBC),
1580  48*24*22,0.,48*24*22) );
1581  for (Int_t iCh = 0; iCh < 48*24*22; iCh++)
1582  SetEMCALChannelTimeRecalibrationFactor(iBC,iCh,0.,kTRUE);
1583  }
1584  }
1585 
1586 
1587 
1588  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1589  fEMCALTimeRecalibrationFactors->Compress();
1590 
1591  // In order to avoid rewriting the same histograms
1592  TH1::AddDirectory(oldStatus);
1593 }
1594 
1598 //____________________________________________________
1600 {
1601  AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1602 
1603  // In order to avoid rewriting the same histograms
1604  Bool_t oldStatus = TH1::AddDirectoryStatus();
1605  TH1::AddDirectory(kFALSE);
1606 
1607  fEMCALBadChannelMap = new TObjArray(22);
1608  //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
1609 
1610  for (int i = 0; i < 22; i++)
1611  fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1612 
1613  fEMCALBadChannelMap->SetOwner(kTRUE);
1614  fEMCALBadChannelMap->Compress();
1615 
1616  // In order to avoid rewriting the same histograms
1617  TH1::AddDirectory(oldStatus);
1618 }
1619 
1623 //___________________________________________________________
1625 {
1626  AliDebug(2,"AliEMCALRecoUtils::InitEMCALL1PhaseInTimeRecalibrationFactors()");
1627 
1628  // In order to avoid rewriting the same histograms
1629  Bool_t oldStatus = TH1::AddDirectoryStatus();
1630  TH1::AddDirectory(kFALSE);
1631 
1633 
1634  fEMCALL1PhaseInTimeRecalibration->Add(new TH1C("h0","EMCALL1phaseForSM", 22, 0, 22));
1635 
1636  for (Int_t i = 0; i < 22; i++) //loop over SMs, default value = 0
1638 
1639  fEMCALL1PhaseInTimeRecalibration->SetOwner(kTRUE);
1641 
1642  // In order to avoid rewriting the same histograms
1643  TH1::AddDirectory(oldStatus);
1644 }
1645 
1655 //____________________________________________________________________________
1656 void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
1657  AliVCluster * cluster,
1658  AliVCaloCells * cells,
1659  Int_t bc)
1660 {
1661  if (!cluster)
1662  {
1663  AliInfo("Cluster pointer null!");
1664  return;
1665  }
1666 
1667  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1668  UShort_t * index = cluster->GetCellsAbsId() ;
1669  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1670  Int_t ncells = cluster->GetNCells();
1671 
1672  // Initialize some used variables
1673  Float_t energy = 0;
1674  Int_t absId =-1;
1675  Int_t icol =-1, irow =-1, imod=1;
1676  Float_t factor = 1, frac = 0;
1677  Int_t absIdMax = -1;
1678  Float_t emax = 0;
1679 
1680  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1681  for (Int_t icell = 0; icell < ncells; icell++)
1682  {
1683  absId = index[icell];
1684  frac = fraction[icell];
1685  if (frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
1686 
1688  {
1689  // Energy
1690  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1691  geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1692  if (fEMCALRecalibrationFactors->GetEntries() <= imod)
1693  continue;
1694  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1695  factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1696 
1697  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1698  imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1699 
1700  }
1701 
1702  energy += cells->GetCellAmplitude(absId)*factor*frac;
1703 
1704  if (emax < cells->GetCellAmplitude(absId)*factor*frac)
1705  {
1706  emax = cells->GetCellAmplitude(absId)*factor*frac;
1707  absIdMax = absId;
1708  }
1709  }
1710 
1711  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
1712 
1713  cluster->SetE(energy);
1714 
1715  // Recalculate time of cluster
1716  Double_t timeorg = cluster->GetTOF();
1717  Bool_t isLowGain = !(cells->GetCellHighGain(absIdMax));//HG = false -> LG = true
1718 
1719  Double_t time = cells->GetCellTime(absIdMax);
1721  RecalibrateCellTime(absIdMax,bc,time,isLowGain);
1722  time-=fConstantTimeShift*1e-9; // only in case of Run1 old simulations
1723 
1724  // Recalibrate time with L1 phase
1726  RecalibrateCellTimeL1Phase(imod, bc, time);
1727 
1728  cluster->SetTOF(time);
1729 
1730  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
1731 }
1732 
1740 //_______________________________________________________________________
1741 void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, Int_t bc)
1742 {
1744  return;
1745 
1746  if (!cells)
1747  {
1748  AliInfo("Cells pointer null!");
1749  return;
1750  }
1751 
1752  Short_t absId =-1;
1753  Bool_t accept = kFALSE;
1754  Float_t ecell = 0;
1755  Double_t tcell = 0;
1756  Double_t ecellin = 0;
1757  Double_t tcellin = 0;
1758  Int_t mclabel = -1;
1759  Double_t efrac = 0;
1760 
1761  Int_t nEMcell = cells->GetNumberOfCells() ;
1762  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1763  {
1764  cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
1765 
1766  accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
1767  if (!accept)
1768  {
1769  ecell = 0;
1770  tcell = -1;
1771  }
1772 
1773  // Set new values
1774  cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
1775  }
1776 
1777  fCellsRecalibrated = kTRUE;
1778 }
1779 
1788 //_______________________________________________________________________________________________________
1789 void AliEMCALRecoUtils::RecalibrateCellTime(Int_t absId, Int_t bc, Double_t & celltime, Bool_t isLGon) const
1790 {
1791  if (!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0) {
1792  if(fLowGain)
1793  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,isLGon)*1.e-9;
1794  else
1795  celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId,kFALSE)*1.e-9;
1796  }
1797 }
1798 
1806 //_______________________________________________________________________________________________________
1808 {
1809  if (!fCellsRecalibrated && IsL1PhaseInTimeRecalibrationOn() && bc >= 0)
1810  {
1811  bc=bc%4;
1812 
1813  Float_t offsetPerSM=0.;
1814  Int_t l1PhaseShift = GetEMCALL1PhaseInTimeRecalibrationForSM(iSM);
1815  Int_t l1Phase=l1PhaseShift & 3; //bit operation
1816 
1817  if(bc >= l1Phase)
1818  offsetPerSM = (bc - l1Phase)*25;
1819  else
1820  offsetPerSM = (bc - l1Phase + 4)*25;
1821 
1822  Int_t l1shiftOffset=l1PhaseShift>>2; //bit operation
1823  l1shiftOffset*=25;
1824 
1825  celltime -= offsetPerSM*1.e-9;
1826  celltime -= l1shiftOffset*1.e-9;
1827  }
1828 }
1829 
1845 //______________________________________________________________________________
1846 void AliEMCALRecoUtils::RecalculateCellLabelsRemoveAddedGenerator( Int_t absID, AliVCluster* clus, AliMCEvent* mc,
1847  Float_t & amp, TArrayI & labeArr, TArrayF & eDepArr) const
1848 {
1849  TString genName ;
1850  Float_t eDepFrac[4];
1851 
1852  Float_t edepTotFrac = 1;
1853  Bool_t found = kFALSE;
1854  Float_t ampOrg = amp;
1855 
1856  //
1857  // Get the energy deposition fraction from cluster.
1858  //
1859  for(Int_t icluscell = 0; icluscell < clus->GetNCells(); icluscell++ )
1860  {
1861  if(absID == clus->GetCellAbsId(icluscell))
1862  {
1863  clus->GetCellMCEdepFractionArray(icluscell,eDepFrac);
1864 
1865  found = kTRUE;
1866 
1867  break;
1868  }
1869  }
1870 
1871  if ( !found )
1872  {
1873  AliWarning(Form("Cell abs ID %d NOT found in cluster",absID));
1874  return;
1875  }
1876 
1877  //
1878  // Check if there is a particle contribution from a given generator name.
1879  // If it is not one of the selected generators,
1880  // remove the constribution from the cell.
1881  //
1882  if ( mc && fNMCGenerToAccept > 0 )
1883  {
1884  //printf("Accept contribution from generator? \n");
1885  for(UInt_t imc = 0; imc < 4; imc++)
1886  {
1887  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc )
1888  {
1889  mc->GetCocktailGenerator(clus->GetLabelAt(imc),genName);
1890 
1891  Bool_t generOK = kFALSE;
1892  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
1893  {
1894  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
1895  }
1896 
1897  if ( !generOK )
1898  {
1899  amp-=ampOrg*eDepFrac[imc];
1900 
1901  edepTotFrac-=eDepFrac[imc];
1902 
1903  eDepFrac[imc] = 0;
1904  }
1905 
1906  } // eDep > 0
1907  } // 4 possible loop
1908 
1909  } // accept at least one kind of generator
1910 
1911  //
1912  // Add MC label and Edep to corresponding array (to be used later in digits)
1913  //
1914  Int_t nLabels = 0;
1915  for(UInt_t imc = 0; imc < 4; imc++)
1916  {
1917  if ( eDepFrac[imc] > 0 && clus->GetNLabels() > imc && edepTotFrac > 0 )
1918  {
1919  nLabels++;
1920 
1921  labeArr.Set(nLabels);
1922  labeArr.AddAt(clus->GetLabelAt(imc), nLabels-1);
1923 
1924  eDepArr.Set(nLabels);
1925  eDepArr.AddAt( (eDepFrac[imc]/edepTotFrac) * amp, nLabels-1);
1926  // use as deposited energy a fraction of the simulated energy (smeared and with noise)
1927  } // edep > 0
1928 
1929  } // mc cell label loop
1930 
1931  //
1932  // If no label found, reject cell
1933  // It can happen to have this case (4 MC labels per cell is not enough for some cases)
1934  // better to remove. To be treated carefully.
1935  //
1936  if ( nLabels == 0 ) amp = 0;
1937 
1938 }
1939 
1947 //______________________________________________________________________________
1948 void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
1949  AliVCaloCells* cells,
1950  AliVCluster* clu)
1951 {
1952  if (!clu)
1953  {
1954  AliInfo("Cluster pointer null!");
1955  return;
1956  }
1957 
1959  else if (fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1960  else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
1961 }
1962 
1971 //_____________________________________________________________________________________________
1973  AliVCaloCells* cells,
1974  AliVCluster* clu)
1975 {
1976  Double_t eCell = 0.;
1977  Float_t fraction = 1.;
1978  Float_t recalFactor = 1.;
1979 
1980  Int_t absId = -1;
1981  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1982  Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
1983  Float_t weight = 0., totalWeight=0.;
1984  Float_t newPos[3] = {-1.,-1.,-1.};
1985  Double_t pLocal[3], pGlobal[3];
1986  Bool_t shared = kFALSE;
1987 
1988  Float_t clEnergy = clu->E(); //Energy already recalibrated previously
1989  if (clEnergy <= 0) return;
1990 
1991  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
1992  Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1993 
1994  //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1995 
1996  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1997  {
1998  absId = clu->GetCellAbsId(iDig);
1999  fraction = clu->GetCellAmplitudeFraction(iDig);
2000  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2001 
2002  if (!fCellsRecalibrated)
2003  {
2004  geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
2005  geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
2006  if (IsRecalibrationOn()) {
2007  recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
2008  }
2009  }
2010 
2011  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2012 
2013  weight = GetCellWeight(eCell,clEnergy);
2014  totalWeight += weight;
2015 
2016  geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
2017  //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
2018  geom->GetGlobal(pLocal,pGlobal,iSupModMax);
2019  //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
2020 
2021  for (int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
2022  }// cell loop
2023 
2024  if (totalWeight>0)
2025  {
2026  for (int i=0; i<3; i++ ) newPos[i] /= totalWeight;
2027  }
2028 
2029  //Float_t pos[]={0,0,0};
2030  //clu->GetPosition(pos);
2031  //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
2032  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
2033 
2034  if (iSupModMax > 1) { //sector 1
2035  newPos[0] +=fMisalTransShift[3];//-=3.093;
2036  newPos[1] +=fMisalTransShift[4];//+=6.82;
2037  newPos[2] +=fMisalTransShift[5];//+=1.635;
2038  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
2039  } else { //sector 0
2040  newPos[0] +=fMisalTransShift[0];//+=1.134;
2041  newPos[1] +=fMisalTransShift[1];//+=8.2;
2042  newPos[2] +=fMisalTransShift[2];//+=1.197;
2043  //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
2044  }
2045  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
2046 
2047  clu->SetPosition(newPos);
2048 }
2049 
2058 //____________________________________________________________________________________________
2060  AliVCaloCells* cells,
2061  AliVCluster* clu)
2062 {
2063  Double_t eCell = 1.;
2064  Float_t fraction = 1.;
2065  Float_t recalFactor = 1.;
2066 
2067  Int_t absId = -1;
2068  Int_t iTower = -1;
2069  Int_t iIphi = -1, iIeta = -1;
2070  Int_t iSupMod = -1, iSupModMax = -1;
2071  Int_t iphi = -1, ieta =-1;
2072  Bool_t shared = kFALSE;
2073 
2074  Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
2075 
2076  if (clEnergy <= 0)
2077  return;
2078  GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
2079  Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
2080 
2081  Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
2082  Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
2083  Int_t startingSM = -1;
2084 
2085  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
2086  {
2087  absId = clu->GetCellAbsId(iDig);
2088  fraction = clu->GetCellAmplitudeFraction(iDig);
2089  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2090 
2091  if (iDig==0) startingSM = iSupMod;
2092  else if (iSupMod != startingSM) areInSameSM = kFALSE;
2093 
2094  eCell = cells->GetCellAmplitude(absId);
2095 
2096  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2097  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
2098 
2099  if (!fCellsRecalibrated)
2100  {
2101  if (IsRecalibrationOn()) {
2102  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2103  }
2104  }
2105 
2106  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2107 
2108  weight = GetCellWeight(eCell,clEnergy);
2109  if (weight < 0) weight = 0;
2110  totalWeight += weight;
2111  weightedCol += ieta*weight;
2112  weightedRow += iphi*weight;
2113 
2114  //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
2115  }// cell loop
2116 
2117  Float_t xyzNew[]={-1.,-1.,-1.};
2118  if (areInSameSM == kTRUE)
2119  {
2120  //printf("In Same SM\n");
2121  weightedCol = weightedCol/totalWeight;
2122  weightedRow = weightedRow/totalWeight;
2123  geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2124  }
2125  else
2126  {
2127  //printf("In Different SM\n");
2128  geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
2129  }
2130 
2131  clu->SetPosition(xyzNew);
2132 }
2133 
2141 //___________________________________________________________________________________________
2143  AliVCaloCells* cells,
2144  AliVCluster * cluster)
2145 {
2146  if (!fRecalDistToBadChannels) return;
2147 
2148  if (!cluster)
2149  {
2150  AliInfo("Cluster pointer null!");
2151  return;
2152  }
2153 
2154  // Get channels map of the supermodule where the cluster is.
2155  Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
2156  Bool_t shared = kFALSE;
2157  GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
2158  TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
2159 
2160  Int_t dRrow, dRcol;
2161  Float_t minDist = 10000.;
2162  Float_t dist = 0.;
2163 
2164  // Loop on tower status map
2165  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2166  {
2167  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2168  {
2169  // Check if tower is bad.
2170  if (hMap->GetBinContent(icol,irow)==0) continue;
2171  //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
2172  // iSupMod,icol, irow, icolM,irowM);
2173 
2174  dRrow=TMath::Abs(irowM-irow);
2175  dRcol=TMath::Abs(icolM-icol);
2176  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2177  if (dist < minDist)
2178  {
2179  //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
2180  minDist = dist;
2181  }
2182  }
2183  }
2184 
2185  // In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
2186  if (shared)
2187  {
2188  TH2D* hMap2 = 0;
2189  Int_t iSupMod2 = -1;
2190 
2191  // The only possible combinations are (0,1), (2,3) ... (8,9)
2192  if (iSupMod%2) iSupMod2 = iSupMod-1;
2193  else iSupMod2 = iSupMod+1;
2194  hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
2195 
2196  // Loop on tower status map of second super module
2197  for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
2198  {
2199  for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
2200  {
2201  // Check if tower is bad.
2202  if (hMap2->GetBinContent(icol,irow)==0) continue;
2203 
2204  //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",
2205  // iSupMod2,icol, irow,iSupMod,icolM,irowM);
2206  dRrow=TMath::Abs(irow-irowM);
2207 
2208  if (iSupMod%2)
2209  dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
2210  else
2211  dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
2212 
2213  dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
2214  if (dist < minDist) minDist = dist;
2215  }
2216  }
2217  }// shared cluster in 2 SuperModules
2218 
2219  AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
2220  cluster->SetDistanceToBadChannel(minDist);
2221 }
2222 
2229 //__________________________________________________________________
2230 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
2231 {
2232  if (!cluster)
2233  {
2234  AliInfo("Cluster pointer null!");
2235  return;
2236  }
2237 
2238  if (cluster->GetM02() != 0)
2239  fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
2240 
2241  Float_t pidlist[AliPID::kSPECIESCN+1];
2242  for (Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
2243 
2244  cluster->SetPID(pidlist);
2245 }
2246 
2267 //___________________________________________________________________________________________________________________
2269  AliVCaloCells* cells, AliVCluster * cluster,
2270  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2271  Float_t & enAfterCuts, Float_t & l0, Float_t & l1,
2272  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2273  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2274 {
2275  if (!cluster)
2276  {
2277  AliInfo("Cluster pointer null!");
2278  return;
2279  }
2280 
2281  Double_t eCell = 0.;
2282  Double_t tCell = 0.;
2283  Float_t fraction = 1.;
2284  Float_t recalFactor = 1.;
2285  Bool_t isLowGain = kFALSE;
2286 
2287  Int_t iSupMod = -1;
2288  Int_t iTower = -1;
2289  Int_t iIphi = -1;
2290  Int_t iIeta = -1;
2291  Int_t iphi = -1;
2292  Int_t ieta = -1;
2293  Double_t etai = -1.;
2294  Double_t phii = -1.;
2295 
2296  Int_t nstat = 0 ;
2297  Float_t wtot = 0.;
2298  Double_t w = 0.;
2299  Double_t etaMean = 0.;
2300  Double_t phiMean = 0.;
2301 
2302  Double_t pGlobal[3];
2303 
2304  // Loop on cells, calculate the cluster energy, in case a cut on cell energy is added,
2305  // or the non linearity correction was applied
2306  // and to check if the cluster is between 2 SM in eta
2307  Int_t iSM0 = -1;
2308  Bool_t shared = kFALSE;
2309  Float_t energy = 0;
2310 
2311  enAfterCuts = 0;
2312  l0 = 0; l1 = 0;
2313  disp = 0; dEta = 0; dPhi = 0;
2314  sEta = 0; sPhi = 0; sEtaPhi = 0;
2315 
2316  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2317  {
2318  // Get from the absid the supermodule, tower and eta/phi numbers
2319 
2320  Int_t absId = cluster->GetCellAbsId(iDigit);
2321 
2322  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2323  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2324 
2325  // Check if there are cells of different SM
2326  if (iDigit == 0 ) iSM0 = iSupMod;
2327  else if (iSupMod!= iSM0) shared = kTRUE;
2328 
2329  // Get the cell energy, if recalibration is on, apply factors
2330  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2331  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2332 
2333  if (IsRecalibrationOn())
2334  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2335 
2336  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2337  tCell = cells->GetCellTime (absId);
2338  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2339 
2340  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2341  tCell*=1e9;
2342  tCell-=fConstantTimeShift;
2343 
2344  if(eCell > 0.05) // at least a minimum 50 MeV cut
2345  energy += eCell;
2346 
2347  if(eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut)
2348  enAfterCuts += eCell;
2349  } // cell loop
2350 
2351 
2352  // Loop on cells to calculate weights and shower shape terms parameters
2353  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2354  {
2355  // Get from the absid the supermodule, tower and eta/phi numbers
2356  Int_t absId = cluster->GetCellAbsId(iDigit);
2357 
2358  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2359  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2360 
2361  //Get the cell energy, if recalibration is on, apply factors
2362  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2363  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2364 
2366  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2367 
2368  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2369  tCell = cells->GetCellTime (absId);
2370  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2371 
2372  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2373  tCell*=1e9;
2374  tCell-=fConstantTimeShift;
2375 
2376  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2377  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2378  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2379 
2380  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2381  {
2382  w = GetCellWeight(eCell, energy);
2383 
2384  // Cell index
2385  if ( fShowerShapeCellLocationType == 0 )
2386  {
2387  etai=(Double_t)ieta;
2388  phii=(Double_t)iphi;
2389  }
2390  // Cell angle location
2391  else if( fShowerShapeCellLocationType == 1 )
2392  {
2393  geom->EtaPhiFromIndex(absId, etai, phii);
2394  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2395  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2396  }
2397  else
2398  {
2399  geom->GetGlobal(absId,pGlobal);
2400 
2401  // Cell x-z location
2402  if( fShowerShapeCellLocationType == 2 )
2403  {
2404  etai = pGlobal[2];
2405  phii = pGlobal[0];
2406  }
2407  // Cell r-z location
2408  else
2409  {
2410  etai = pGlobal[2];
2411  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2412  }
2413  }
2414 
2415  if (w > 0.0)
2416  {
2417  wtot += w ;
2418  nstat++;
2419 
2420  // Shower shape
2421  sEta += w * etai * etai ;
2422  etaMean += w * etai ;
2423  sPhi += w * phii * phii ;
2424  phiMean += w * phii ;
2425  sEtaPhi += w * etai * phii ;
2426  }
2427  }
2428  else if(eCell > 0.05)
2429  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2430  } // cell loop
2431 
2432  //printf("sEta %f sPhi %f etaMean %f phiMean %f sEtaPhi %f wtot %f\n",sEta,sPhi,etaMean,phiMean,sEtaPhi, wtot);
2433 
2434  // Normalize to the weight
2435  if (wtot > 0)
2436  {
2437  etaMean /= wtot ;
2438  phiMean /= wtot ;
2439  }
2440  else
2441  AliDebug(2,Form("Wrong weight %f\n", wtot));
2442 
2443  // Loop on cells to calculate dispersion
2444  for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2445  {
2446  // Get from the absid the supermodule, tower and eta/phi numbers
2447  Int_t absId = cluster->GetCellAbsId(iDigit);
2448 
2449  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
2450  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2451 
2452  //Get the cell energy, if recalibration is on, apply factors
2453  fraction = cluster->GetCellAmplitudeFraction(iDigit);
2454  if (fraction < 1e-4) fraction = 1.; // in case unfolding is off
2455  if (IsRecalibrationOn())
2456  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2457 
2458  eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
2459  tCell = cells->GetCellTime (absId);
2460  isLowGain = !(cells->GetCellHighGain(absId));//HG = false -> LG = true
2461 
2462  RecalibrateCellTime(absId, bc, tCell,isLowGain);
2463  tCell*=1e9;
2464  tCell-=fConstantTimeShift;
2465 
2466  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2467  // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2468  if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2469 
2470  if ( energy > 0 && eCell > cellEcut && TMath::Abs(tCell) < cellTimeCut )
2471  {
2472  w = GetCellWeight(eCell,cluster->E());
2473 
2474  // Cell index
2475  if ( fShowerShapeCellLocationType == 0 )
2476  {
2477  etai=(Double_t)ieta;
2478  phii=(Double_t)iphi;
2479  }
2480  // Cell angle location
2481  else if( fShowerShapeCellLocationType == 1 )
2482  {
2483  geom->EtaPhiFromIndex(absId, etai, phii);
2484  etai *= TMath::RadToDeg(); // change units to degrees instead of radians
2485  phii *= TMath::RadToDeg(); // change units to degrees instead of radians
2486  }
2487  else
2488  {
2489  geom->GetGlobal(absId,pGlobal);
2490 
2491  // Cell x-z location
2492  if( fShowerShapeCellLocationType == 2 )
2493  {
2494  etai = pGlobal[2];
2495  phii = pGlobal[0];
2496  }
2497  // Cell r-z location
2498  else
2499  {
2500  etai = pGlobal[2];
2501  phii = TMath::Sqrt(pGlobal[0]*pGlobal[0]+pGlobal[1]*pGlobal[1]);
2502  }
2503  }
2504 
2505  if (w > 0.0)
2506  {
2507  disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
2508  dEta += w * (etai-etaMean)*(etai-etaMean) ;
2509  dPhi += w * (phii-phiMean)*(phii-phiMean) ;
2510  }
2511  }
2512  else
2513  AliDebug(2,Form("Wrong energy in cell %f and/or cluster %f\n", eCell, cluster->E()));
2514  }// cell loop
2515 
2516  // Normalize to the weigth and set shower shape parameters
2517  if (wtot > 0 && nstat > 1)
2518  {
2519  disp /= wtot ;
2520  dEta /= wtot ;
2521  dPhi /= wtot ;
2522  sEta /= wtot ;
2523  sPhi /= wtot ;
2524  sEtaPhi /= wtot ;
2525 
2526  sEta -= etaMean * etaMean ;
2527  sPhi -= phiMean * phiMean ;
2528  sEtaPhi -= etaMean * phiMean ;
2529 
2530  l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2531  l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2532 
2533  //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);
2534 
2535  }
2536  else
2537  {
2538  l0 = 0. ;
2539  l1 = 0. ;
2540  dEta = 0. ; dPhi = 0. ; disp = 0. ;
2541  sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
2542  }
2543 }
2544 
2563 //___________________________________________________________________________________________________________________
2565  AliVCaloCells* cells, AliVCluster * cluster,
2566  Float_t & l0, Float_t & l1,
2567  Float_t & disp, Float_t & dEta, Float_t & dPhi,
2568  Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
2569 {
2570  Float_t newEnergy = 0;
2571  Float_t cellEmin = 0.05; // 50 MeV
2572  Float_t cellTimeWindow = 1000; // open cut
2573  Int_t bc = 0;
2575  cellEmin, cellTimeWindow, bc,
2576  newEnergy, l0, l1, disp,
2577  dEta, dPhi, sEta, sPhi, sEtaPhi);
2578 }
2579 
2588 //____________________________________________________________________________________________
2590  AliVCaloCells* cells,
2591  AliVCluster * cluster)
2592 {
2593  Float_t l0 = 0., l1 = 0.;
2594  Float_t disp = 0., dEta = 0., dPhi = 0.;
2595  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2596 
2597  AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
2598  dEta, dPhi, sEta, sPhi, sEtaPhi);
2599 
2600  cluster->SetM02(l0);
2601  cluster->SetM20(l1);
2602  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2603 
2604 }
2605 
2618 //____________________________________________________________________________________________
2620  AliVCaloCells* cells, AliVCluster * cluster,
2621  Float_t cellEcut, Float_t cellTimeCut, Int_t bc,
2622  Float_t & enAfterCuts)
2623 {
2624  Float_t l0 = 0., l1 = 0.;
2625  Float_t disp = 0., dEta = 0., dPhi = 0.;
2626  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2627 
2629  cellEcut, cellTimeCut, bc,
2630  enAfterCuts, l0, l1, disp,
2631  dEta, dPhi, sEta, sPhi, sEtaPhi);
2632 
2633  cluster->SetM02(l0);
2634  cluster->SetM20(l1);
2635  if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
2636 
2637 }
2638 
2652 //____________________________________________________________________________
2653 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
2654  TObjArray * clusterArr,
2655  const AliEMCALGeometry *geom,
2656  AliMCEvent * mc)
2657 {
2658  fMatchedTrackIndex ->Reset();
2659  fMatchedClusterIndex->Reset();
2660  fResidualPhi->Reset();
2661  fResidualEta->Reset();
2662 
2663  fMatchedTrackIndex ->Set(1000);
2664  fMatchedClusterIndex->Set(1000);
2665  fResidualPhi->Set(1000);
2666  fResidualEta->Set(1000);
2667 
2668  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
2669  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
2670 
2671  // Init the magnetic field if not already on
2672  if (!TGeoGlobalMagField::Instance()->GetField())
2673  {
2674  if (!event->InitMagneticField())
2675  {
2676  AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
2677  }
2678  } // Init mag field
2679 
2680  if (esdevent)
2681  {
2682  UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
2683  UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
2684  Bool_t desc1 = (mask1 >> 3) & 0x1;
2685  Bool_t desc2 = (mask2 >> 3) & 0x1;
2686  if (desc1==0 || desc2==0) {
2687  // AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
2688  // mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
2689  // mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
2690  fITSTrackSA=kTRUE;
2691  }
2692  }
2693 
2694  TObjArray *clusterArray = 0x0;
2695  if (!clusterArr)
2696  {
2697  clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
2698  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2699  {
2700  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2701  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells()))
2702  continue;
2703  clusterArray->AddAt(cluster,icl);
2704  }
2705  }
2706 
2707  Int_t matched=0;
2708  Double_t cv[21];
2709  TString genName;
2710  for (Int_t i=0; i<21;i++) cv[i]=0;
2711  for (Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
2712  {
2713  AliExternalTrackParam *trackParam = 0;
2714  Int_t mcLabel = -1;
2715  // If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
2716  AliESDtrack *esdTrack = 0;
2717  AliAODTrack *aodTrack = 0;
2718  if (esdevent)
2719  {
2720  esdTrack = esdevent->GetTrack(itr);
2721  if (!esdTrack) continue;
2722  if (!IsAccepted(esdTrack)) continue;
2723  if (esdTrack->Pt()<fCutMinTrackPt) continue;
2724 
2725  if ( TMath::Abs(esdTrack->Eta()) > 0.9 ) continue;
2726 
2727  // Save some time and memory in case of no DCal present
2728  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2729  {
2730  Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
2731  if ( phi <= 10 || phi >= 250 ) continue;
2732  }
2733 
2734  if (!fITSTrackSA)
2735  trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); // if TPC Available
2736  else
2737  trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
2738 
2739  mcLabel = TMath::Abs(esdTrack->GetLabel());
2740  }
2741 
2742  // If the input event is AOD, the starting point for extrapolation is at vertex
2743  // AOD tracks are selected according to its filterbit.
2744  else if (aodevent)
2745  {
2746  aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(itr));
2747 
2748  if (!aodTrack) AliFatal("Not a standard AOD");
2749 
2750  if (!aodTrack) continue;
2751 
2752  if (fAODTPCOnlyTracks)
2753  { // Match with TPC only tracks, default from May 2013, before filter bit 32
2754  //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
2755  if (!aodTrack->IsTPCConstrained()) continue ;
2756  }
2757  else if (fAODHybridTracks)
2758  { // Match with hybrid tracks
2759  //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
2760  if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
2761  } else
2762  { // Match with tracks on a mask
2763  //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask));
2764  if (!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks
2765  }
2766 
2767  if (aodTrack->Pt() < fCutMinTrackPt) continue;
2768 
2769  if ( TMath::Abs(aodTrack->Eta()) > 0.9 ) continue;
2770 
2771  // Save some time and memory in case of no DCal present
2772  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2773  {
2774  Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
2775  if ( phi <= 10 || phi >= 250 ) continue;
2776  }
2777 
2778  Double_t pos[3],mom[3];
2779  aodTrack->GetXYZ(pos);
2780  aodTrack->GetPxPyPz(mom);
2781  AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",
2782  itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
2783 
2784  trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
2785 
2786  mcLabel = TMath::Abs(aodTrack->GetLabel());
2787  }
2788 
2789  //Return if the input data is not "AOD" or "ESD"
2790  else
2791  {
2792  AliWarning("Wrong input data type! Should be \"AOD\" or \"ESD\" ");
2793  if (clusterArray)
2794  {
2795  clusterArray->Clear();
2796  delete clusterArray;
2797  }
2798  return;
2799  }
2800 
2801  if (!trackParam) continue;
2802 
2803  //
2804  // Check if track comes from a particular MC generator, do not include it
2805  // if it is not a selected one
2806  //
2807  if( mc && fMCGenerToAcceptForTrack && fNMCGenerToAccept > 0 )
2808  {
2809  mc->GetCocktailGenerator(mcLabel,genName);
2810 
2811  Bool_t generOK = kFALSE;
2812  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
2813  {
2814  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
2815  }
2816 
2817  if ( !generOK ) continue;
2818  }
2819 
2820  // Extrapolate the track to EMCal surface, see AliEMCALRecoUtilsBase
2821  AliExternalTrackParam emcalParam(*trackParam);
2822  Float_t eta, phi, pt;
2823  if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2824  {
2825  if (aodevent && trackParam) delete trackParam;
2826  if (fITSTrackSA && trackParam) delete trackParam;
2827  continue;
2828  }
2829 
2830  if ( TMath::Abs(eta) > 0.75 )
2831  {
2832  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2833  continue;
2834  }
2835 
2836  // Save some time and memory in case of no DCal present
2837  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2838  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2839  {
2840  if ( trackParam && (aodevent || fITSTrackSA) ) delete trackParam;
2841  continue;
2842  }
2843 
2844  //Find matched clusters
2845  Int_t index = -1;
2846  Float_t dEta = -999, dPhi = -999;
2847  if (!clusterArr)
2848  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
2849  else
2850  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr , dEta, dPhi);
2851 
2852 
2853  if (index>-1)
2854  {
2855  fMatchedTrackIndex ->AddAt(itr,matched);
2856  fMatchedClusterIndex ->AddAt(index,matched);
2857  fResidualEta ->AddAt(dEta,matched);
2858  fResidualPhi ->AddAt(dPhi,matched);
2859  matched++;
2860  }
2861 
2862  if (aodevent && trackParam) delete trackParam;
2863  if (fITSTrackSA && trackParam) delete trackParam;
2864  }//track loop
2865 
2866  if (clusterArray)
2867  {
2868  clusterArray->Clear();
2869  delete clusterArray;
2870  }
2871 
2872  AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
2873 
2874  fMatchedTrackIndex ->Set(matched);
2875  fMatchedClusterIndex ->Set(matched);
2876  fResidualPhi ->Set(matched);
2877  fResidualEta ->Set(matched);
2878 }
2879 
2891 //________________________________________________________________________________
2893  const AliVEvent *event,
2894  const AliEMCALGeometry *geom,
2895  Float_t &dEta, Float_t &dPhi)
2896 {
2897  Int_t index = -1;
2898 
2899  if ( TMath::Abs(track->Eta()) > 0.9 ) return index;
2900 
2901  // Save some time and memory in case of no DCal present
2902  if( geom->GetNumberOfSuperModules() < 13 ) // Run1 10 (12, 2 not active but present)
2903  {
2904  Double_t phiV = track->Phi()*TMath::RadToDeg();
2905  if ( phiV <= 10 || phiV >= 250 ) return index;
2906  }
2907 
2908  AliExternalTrackParam *trackParam = 0;
2909  if (!fITSTrackSA)
2910  trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam()); // If TPC
2911  else
2912  trackParam = new AliExternalTrackParam(*track);
2913 
2914  if (!trackParam) return index;
2915  AliExternalTrackParam emcalParam(*trackParam);
2916 
2917  Float_t eta, phi, pt;
2918  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt))
2919  {
2920  if (fITSTrackSA) delete trackParam;
2921  return index;
2922  }
2923 
2924  if ( TMath::Abs(eta) > 0.75 )
2925  {
2926  if (fITSTrackSA) delete trackParam;
2927  return index;
2928  }
2929 
2930  // Save some time and memory in case of no DCal present
2931  if ( geom->GetNumberOfSuperModules() < 13 && // Run1 10 (12, 2 not active but present)
2932  ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
2933  {
2934  if (fITSTrackSA) delete trackParam;
2935  return index;
2936  }
2937 
2938  TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
2939 
2940  for (Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2941  {
2942  AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2943  if (!IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
2944  clusterArr->AddAt(cluster,icl);
2945  }
2946 
2947  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
2948  clusterArr->Clear();
2949  delete clusterArr;
2950  if (fITSTrackSA) delete trackParam;
2951 
2952  return index;
2953 }
2954 
2965 //_______________________________________________________________________________________________
2966 Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
2967  AliExternalTrackParam *trkParam,
2968  const TObjArray * clusterArr,
2969  Float_t &dEta, Float_t &dPhi)
2970 {
2971  dEta=-999, dPhi=-999;
2972  Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
2973  Int_t index = -1;
2974  Float_t tmpEta=-999, tmpPhi=-999;
2975 
2976  Double_t exPos[3] = {0.,0.,0.};
2977  if (!emcalParam->GetXYZ(exPos)) return index;
2978 
2979  Float_t clsPos[3] = {0.,0.,0.};
2980  for (Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
2981  {
2982  AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
2983 
2984  if (!cluster || !cluster->IsEMCAL()) continue;
2985 
2986  cluster->GetPosition(clsPos);
2987 
2988  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));
2989  if (dR > fClusterWindow) continue;
2990 
2991  AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
2992 
2993  if (!AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
2994 
2995  if (fCutEtaPhiSum)
2996  {
2997  Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
2998  if (tmpR<dRMax)
2999  {
3000  dRMax=tmpR;
3001  dEtaMax=tmpEta;
3002  dPhiMax=tmpPhi;
3003  index=icl;
3004  }
3005  }
3006  else if (fCutEtaPhiSeparate)
3007  {
3008  if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
3009  {
3010  dEtaMax = tmpEta;
3011  dPhiMax = tmpPhi;
3012  index=icl;
3013  }
3014  }
3015  else
3016  {
3017  AliError("Please specify your cut criteria");
3018  AliError("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()");
3019  AliError("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()");
3020  return index;
3021  }
3022  }
3023 
3024  dEta=dEtaMax;
3025  dPhi=dPhiMax;
3026 
3027  return index;
3028 }
3029 
3042 //---------------------------------------------------------------------------------
3043 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
3044  const AliVCluster *cluster,
3045  Float_t &tmpEta,
3046  Float_t &tmpPhi)
3047 {
3048  return AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
3049 }
3050 
3060 //_______________________________________________________________________
3062  Float_t &dEta, Float_t &dPhi)
3063 {
3064  if (FindMatchedPosForCluster(clsIndex) >= 999)
3065  {
3066  AliDebug(2,"No matched tracks found!\n");
3067  dEta=999.;
3068  dPhi=999.;
3069  return;
3070  }
3071 
3072  dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
3073  dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
3074 }
3075 
3085 //______________________________________________________________________________________________
3087 {
3088  if (FindMatchedPosForTrack(trkIndex) >= 999)
3089  {
3090  AliDebug(2,"No matched cluster found!\n");
3091  dEta=999.;
3092  dPhi=999.;
3093  return;
3094  }
3095 
3096  dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
3097  dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
3098 }
3099 
3106 //__________________________________________________________
3108 {
3109  if (IsClusterMatched(clsIndex))
3110  return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
3111  else
3112  return -1;
3113 }
3114 
3121 //__________________________________________________________
3123 {
3124  if (IsTrackMatched(trkIndex))
3125  return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
3126  else
3127  return -1;
3128 }
3129 
3137 //______________________________________________________________
3139 {
3140  if (FindMatchedPosForCluster(clsIndex) < 999)
3141  return kTRUE;
3142  else
3143  return kFALSE;
3144 }
3145 
3153 //____________________________________________________________
3155 {
3156  if (FindMatchedPosForTrack(trkIndex) < 999)
3157  return kTRUE;
3158  else
3159  return kFALSE;
3160 }
3161 
3169 //______________________________________________________________________
3171 {
3172  Float_t tmpR = fCutR;
3173  UInt_t pos = 999;
3174 
3175  for (Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
3176  {
3177  if (fMatchedClusterIndex->At(i)==clsIndex)
3178  {
3179  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3180  if (r<tmpR)
3181  {
3182  pos=i;
3183  tmpR=r;
3184  AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3185  fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3186  }
3187  }
3188  }
3189  return pos;
3190 }
3191 
3199 //____________________________________________________________________
3201 {
3202  Float_t tmpR = fCutR;
3203  UInt_t pos = 999;
3204 
3205  for (Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
3206  {
3207  if (fMatchedTrackIndex->At(i)==trkIndex)
3208  {
3209  Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
3210  if (r<tmpR)
3211  {
3212  pos=i;
3213  tmpR=r;
3214  AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
3215  fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
3216  }
3217  }
3218  }
3219  return pos;
3220 }
3221 
3236 //__________________________________________________________________________
3238  const AliEMCALGeometry *geom,
3239  AliVCaloCells* cells, Int_t bc)
3240 {
3241  Bool_t isGood=kTRUE;
3242 
3243  if (!cluster || !cluster->IsEMCAL()) return kFALSE;
3244  if (ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
3245  if (!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
3246  if (IsExoticCluster(cluster, cells,bc)) return kFALSE;
3247 
3248  return isGood;
3249 }
3250 
3261 //__________________________________________________________
3263 {
3264  UInt_t status = esdTrack->GetStatus();
3265 
3266  Int_t nClustersITS = esdTrack->GetITSclusters(0);
3267  Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
3268 
3269  Float_t chi2PerClusterITS = -1;
3270  Float_t chi2PerClusterTPC = -1;
3271  if (nClustersITS!=0)
3272  chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
3273  if (nClustersTPC!=0)
3274  chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
3275 
3276  //
3277  // DCA cuts
3278  // Only to be used for primary
3279  //
3280  if ( fTrackCutsType == kGlobalCut )
3281  {
3282  Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01);
3283  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
3284 
3285  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3286 
3287  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3288  }
3289  else if( fTrackCutsType == kGlobalCut2011 )
3290  {
3291  Float_t maxDCAToVertexXYPtDep = 0.0105 + 0.0350/TMath::Power(esdTrack->Pt(),1.1);
3292  // This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2011()
3293 
3294  //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
3295 
3296  SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
3297  }
3298 
3299  Float_t b[2];
3300  Float_t bCov[3];
3301  esdTrack->GetImpactParameters(b,bCov);
3302  if (bCov[0]<=0 || bCov[2]<=0)
3303  {
3304  AliDebug(1, "Estimated b resolution lower or equal zero!");
3305  bCov[0]=0; bCov[2]=0;
3306  }
3307 
3308  Float_t dcaToVertexXY = b[0];
3309  Float_t dcaToVertexZ = b[1];
3310  Float_t dcaToVertex = -1;
3311 
3312  if (fCutDCAToVertex2D)
3313  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY / fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY +
3314  dcaToVertexZ*dcaToVertexZ / fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ );
3315  else
3316  dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
3317 
3318  // cut the track?
3319 
3320  Bool_t cuts[kNCuts];
3321  for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
3322 
3323  // track quality cuts
3324  if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
3325  cuts[0]=kTRUE;
3326  if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
3327  cuts[1]=kTRUE;
3328  if (nClustersTPC<fCutMinNClusterTPC)
3329  cuts[2]=kTRUE;
3330  if (nClustersITS<fCutMinNClusterITS)
3331  cuts[3]=kTRUE;
3332  if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
3333  cuts[4]=kTRUE;
3334  if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
3335  cuts[5]=kTRUE;
3336  if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
3337  cuts[6]=kTRUE;
3338  if (fCutDCAToVertex2D && dcaToVertex > 1)
3339  cuts[7] = kTRUE;
3340  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
3341  cuts[8] = kTRUE;
3342  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
3343  cuts[9] = kTRUE;
3344 
3346  {
3347  //Require at least one SPD point + anything else in ITS
3348  if ( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
3349  cuts[10] = kTRUE;
3350  }
3351 
3352  // ITS
3354  {
3355  if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin))
3356  {
3357  // TPC tracks
3358  cuts[11] = kTRUE;
3359  }
3360  else
3361  {
3362  // ITS standalone tracks
3364  {
3365  if (status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
3366  }
3367  else if (fCutRequireITSpureSA)
3368  {
3369  if (!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
3370  }
3371  }
3372  }
3373 
3374  Bool_t cut=kFALSE;
3375  for (Int_t i=0; i<kNCuts; i++)
3376  if (cuts[i]) { cut = kTRUE ; }
3377 
3378  // cut the track
3379  if (cut)
3380  return kFALSE;
3381  else
3382  return kTRUE;
3383 }
3384 
3390 //_____________________________________
3392 {
3393  switch (fTrackCutsType)
3394  {
3395  case kTPCOnlyCut:
3396  {
3397  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()"));
3398  //TPC
3399  SetMinNClustersTPC(70);
3401  SetAcceptKinkDaughters(kFALSE);
3402  SetRequireTPCRefit(kFALSE);
3403 
3404  //ITS
3405  SetRequireITSRefit(kFALSE);
3406  SetMaxDCAToVertexZ(3.2);
3407  SetMaxDCAToVertexXY(2.4);
3408  SetDCAToVertex2D(kTRUE);
3409 
3410  break;
3411  }
3412 
3413  case kGlobalCut:
3414  {
3415  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE)"));
3416  //TPC
3417  SetMinNClustersTPC(70);
3419  SetAcceptKinkDaughters(kFALSE);
3420  SetRequireTPCRefit(kTRUE);
3421 
3422  //ITS
3423  SetRequireITSRefit(kTRUE);
3424  SetMaxDCAToVertexZ(2);
3426  SetDCAToVertex2D(kFALSE);
3427 
3428  break;
3429  }
3430 
3431  case kLooseCut:
3432  {
3433  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
3434  SetMinNClustersTPC(50);
3435  SetAcceptKinkDaughters(kTRUE);
3436 
3437  break;
3438  }
3439 
3440  case kITSStandAlone:
3441  {
3442  AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
3443  SetRequireITSRefit(kTRUE);
3444  SetRequireITSStandAlone(kTRUE);
3445  SetITSTrackSA(kTRUE);
3446  break;
3447  }
3448 
3449  case kGlobalCut2011:
3450  {
3451  AliInfo(Form("Track cuts for matching: AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE)"));
3452  //TPC
3453  SetMinNClustersTPC(50);
3455  SetAcceptKinkDaughters(kFALSE);
3456  SetRequireTPCRefit(kTRUE);
3457 
3458  //ITS
3459  SetRequireITSRefit(kTRUE);
3460  SetMaxDCAToVertexZ(2);
3462  SetDCAToVertex2D(kFALSE);
3463 
3464  break;
3465  }
3466 
3467  case kLooseCutWithITSrefit:
3468  {
3469  AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut plus ITSrefit"));
3470  SetMinNClustersTPC(50);
3471  SetAcceptKinkDaughters(kTRUE);
3472  SetRequireITSRefit(kTRUE);
3473 
3474  break;
3475  }
3476  }
3477 }
3478 
3484 //________________________________________________________________________
3486 {
3487  Int_t nTracks = event->GetNumberOfTracks();
3488  for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
3489  {
3490  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
3491  if (!track)
3492  {
3493  AliWarning(Form("Could not receive track %d", iTrack));
3494  continue;
3495  }
3496 
3497  Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
3498  track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
3499 
3500  /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
3501  AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
3502  if (esdtrack)
3503  {
3504  if (matchClusIndex != -1)
3505  esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
3506  else
3507  esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3508  }
3509  else
3510  {
3511  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
3512  if (matchClusIndex != -1)
3513  aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
3514  else
3515  aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
3516  }
3517  }
3518  AliDebug(2,"Track matched to closest cluster");
3519 }
3520 
3527 //_________________________________________________________________________
3529 {
3530  for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
3531  {
3532  AliVCluster *cluster = event->GetCaloCluster(iClus);
3533  if (!cluster->IsEMCAL())
3534  continue;
3535 
3536  //
3537  // Remove old matches in cluster
3538  //
3539  if(cluster->GetNTracksMatched() > 0)
3540  {
3541  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
3542  {
3543  TArrayI arrayTrackMatched(0);
3544  ((AliESDCaloCluster*)cluster)->AddTracksMatched(arrayTrackMatched);
3545  }
3546  else
3547  {
3548  for(Int_t iTrack = 0; iTrack < cluster->GetNTracksMatched(); iTrack++)
3549  {
3550  ((AliAODCaloCluster*)cluster)->RemoveTrackMatched((TObject*)((AliAODCaloCluster*)cluster)->GetTrackMatched(iTrack));
3551  }
3552  }
3553  }
3554 
3555  //
3556  // Find new matches and put them in the cluster
3557  //
3558  Int_t nTracks = event->GetNumberOfTracks();
3559  TArrayI arrayTrackMatched(nTracks);
3560 
3561  // Get the closest track matched to the cluster
3562  Int_t nMatched = 0;
3563  Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
3564  if (matchTrackIndex != -1)
3565  {
3566  arrayTrackMatched[nMatched] = matchTrackIndex;
3567  nMatched++;
3568  }
3569 
3570  // Get all other tracks matched to the cluster
3571  for (Int_t iTrk=0; iTrk<nTracks; ++iTrk)
3572  {
3573  AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
3574 
3575  if( !track ) continue;
3576 
3577  if ( iTrk == matchTrackIndex ) continue;
3578 
3579  if ( track->GetEMCALcluster() == iClus )
3580  {
3581  arrayTrackMatched[nMatched] = iTrk;
3582  ++nMatched;
3583  }
3584  }
3585 
3586  AliDebug(2,Form("cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]));
3587 
3588  arrayTrackMatched.Set(nMatched);
3589  AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
3590  if (esdcluster)
3591  esdcluster->AddTracksMatched(arrayTrackMatched);
3592  else if ( nMatched > 0 )
3593  {
3594  AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
3595  if (aodcluster)
3596  {
3597  aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
3598  //AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(event->GetTrack(arrayTrackMatched.At(0)));
3599  //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());
3600  //printf("With specs: pt %.4f, eta %.4f, phi %.4f\n",aodtrack->Pt(),aodtrack->Eta(), aodtrack->Phi());
3601  }
3602  }
3603 
3604  Float_t eta= -999, phi = -999;
3605  if (matchTrackIndex != -1)
3606  GetMatchedResiduals(iClus, eta, phi);
3607 
3608  cluster->SetTrackDistance(phi, eta);
3609  }
3610 
3611  AliDebug(2,"Cluster matched to tracks");
3612 }
3613 
3616  else {
3617  fEMCALRecalibrationFactors = new TObjArray(map->GetEntries());
3618  fEMCALRecalibrationFactors->SetOwner(true);
3619  }
3620  if(!fEMCALRecalibrationFactors->IsOwner()){
3621  // Must claim ownership since the new objects are owend by this instance
3622  fEMCALRecalibrationFactors->SetOwner(kTRUE);
3623  }
3624  for(int i = 0; i < map->GetEntries(); i++){
3625  TH2F *hist = dynamic_cast<TH2F *>(map->At(i));
3626  if(!hist) continue;
3627  this->SetEMCALChannelRecalibrationFactors(i, hist);
3628  }
3629 }
3630 
3634  fEMCALRecalibrationFactors->SetOwner(true);
3635  }
3636  if(fEMCALRecalibrationFactors->GetEntries() <= iSM) fEMCALRecalibrationFactors->Expand(iSM+1);
3637  if(fEMCALRecalibrationFactors->At(iSM)) fEMCALRecalibrationFactors->RemoveAt(iSM);
3638  TH2F *clone = new TH2F(*h);
3639  clone->SetDirectory(NULL);
3640  fEMCALRecalibrationFactors->AddAt(clone,iSM);
3641 }
3642 
3645  else {
3646  fEMCALBadChannelMap = new TObjArray(map->GetEntries());
3647  fEMCALBadChannelMap->SetOwner(true);
3648  }
3649  if(!fEMCALBadChannelMap->IsOwner()) {
3650  // Must claim ownership since the new objects are owend by this instance
3651  fEMCALBadChannelMap->SetOwner(true);
3652  }
3653  for(int i = 0; i < map->GetEntries(); i++){
3654  TH2I *hist = dynamic_cast<TH2I *>(map->At(i));
3655  if(!hist) continue;
3656  this->SetEMCALChannelStatusMap(i, hist);
3657  }
3658 }
3659 
3661  if(!fEMCALBadChannelMap){
3662  fEMCALBadChannelMap = new TObjArray(iSM);
3663  fEMCALBadChannelMap->SetOwner(true);
3664  }
3665  if(fEMCALBadChannelMap->GetEntries() <= iSM) fEMCALBadChannelMap->Expand(iSM+1);
3666  if(fEMCALBadChannelMap->At(iSM)) fEMCALBadChannelMap->RemoveAt(iSM);
3667  TH2I *clone = new TH2I(*h);
3668  clone->SetDirectory(NULL);
3669  fEMCALBadChannelMap->AddAt(clone,iSM);
3670 }
3671 
3674  else {
3675  fEMCALTimeRecalibrationFactors = new TObjArray(map->GetEntries());
3676  fEMCALTimeRecalibrationFactors->SetOwner(true);
3677  }
3678  if(!fEMCALTimeRecalibrationFactors->IsOwner()) {
3679  // Must claim ownership since the new objects are owend by this instance
3680  fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
3681  }
3682  for(int i = 0; i < map->GetEntries(); i++){
3683  TH1F *hist = dynamic_cast<TH1F *>(map->At(i));
3684  if(!hist) continue;
3686  }
3687 }
3688 
3692  fEMCALTimeRecalibrationFactors->SetOwner(true);
3693  }
3694  if(fEMCALTimeRecalibrationFactors->GetEntries() <= bc) fEMCALTimeRecalibrationFactors->Expand(bc+1);
3696  TH1F *clone = new TH1F(*h);
3697  clone->SetDirectory(NULL);
3698  fEMCALTimeRecalibrationFactors->AddAt(clone,bc);
3699 }
3700 
3703  else {
3704  fEMCALL1PhaseInTimeRecalibration = new TObjArray(map->GetEntries());
3705  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3706  }
3707  if(!fEMCALL1PhaseInTimeRecalibration->IsOwner()) {
3708  // Must claim ownership since the new objects are owend by this instance
3709  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3710  }
3711  for(int i = 0; i < map->GetEntries(); i++) {
3712  TH1C *hist = dynamic_cast<TH1C *>(map->At(i));
3713  if(!hist) continue;
3715  }
3716 }
3717 
3721  fEMCALL1PhaseInTimeRecalibration->SetOwner(true);
3722  }
3724  TH1C *clone = new TH1C(*h);
3725  clone->SetDirectory(NULL);
3726  fEMCALL1PhaseInTimeRecalibration->AddAt(clone,0);
3727 }
3728 
3732 //___________________________________________________
3734 {
3735  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3736  printf("AliEMCALRecoUtils Settings: \n");
3737  printf("\tMisalignment shifts\n");
3738  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,
3740  fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
3741  printf("\tNon linearity function %d, parameters:\n", fNonLinearityFunction);
3742  if (fNonLinearityFunction != 3) // print only if not kNoCorrection
3743  for (Int_t i=0; i<10; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
3744 
3745  printf("\tPosition Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
3746 
3747  printf("\tMatching criteria: ");
3748  if (fCutEtaPhiSum) {
3749  printf("\t\tsqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
3750  } else if (fCutEtaPhiSeparate) {
3751  printf("\t\tdEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
3752  } else {
3753  printf("\t\tError\n");
3754  printf("\t\tplease specify your cut criteria\n");
3755  printf("\t\tTo cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
3756  printf("\t\tTo cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
3757  }
3758 
3759  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);
3760  printf("\tCluster selection window: dR < %2.0f\n",fClusterWindow);
3761 
3762  printf("\tTrack cuts: \n");
3763  printf("\t\tMinimum track pT: %1.2f\n",fCutMinTrackPt);
3764  printf("\t\tAOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask);
3765  printf("\t\tTPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
3766  printf("\t\tAcceptKinks = %d\n",fCutAcceptKinkDaughters);
3767  printf("\t\tMinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
3768  printf("\t\tMaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
3769  printf("\t\tDCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
3770  printf("-------------------------------------------------------------------------------------------------------------------------------------- \n");
3771 }
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)