AliPhysics  master (3d17d9d)
AliEmcalCorrectionClusterizer.cxx
Go to the documentation of this file.
1 // AliEmcalCorrectionClusterizer
2 //
3 /**************************************************************************
4  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5  * *
6  * Author: The ALICE Off-line Project. *
7  * Contributors are mentioned in the code where appropriate. *
8  * *
9  * Permission to use, copy, modify and distribute this software and its *
10  * documentation strictly for non-commercial purposes is hereby granted *
11  * without fee, provided that the above copyright notice appears in all *
12  * copies and that both the copyright notice and this permission notice *
13  * appear in the supporting documentation. The authors make no claims *
14  * about the suitability of this software for any purpose. It is *
15  * provided "as is" without express or implied warranty. *
16  **************************************************************************/
17 
18 // --- Root ---
19 #include <TObjArray.h>
20 #include <TArrayI.h>
21 #include <TStopwatch.h>
22 
23 // --- AliRoot ---
24 #include "AliCDBEntry.h"
25 #include "AliCDBManager.h"
26 #include "AliCaloCalibPedestal.h"
27 #include "AliEMCALAfterBurnerUF.h"
28 #include "AliEMCALCalibData.h"
29 #include "AliEMCALClusterizerNxN.h"
30 #include "AliEMCALClusterizerv1.h"
31 #include "AliEMCALClusterizerv2.h"
32 #include "AliEMCALClusterizerv3.h"
33 #include "AliEMCALClusterizerFixedWindow.h"
34 #include "AliEMCALDigit.h"
35 #include "AliEMCALGeometry.h"
36 #include "AliEMCALRecParam.h"
37 #include "AliEMCALRecPoint.h"
38 #include "AliInputEventHandler.h"
39 #include "AliClusterContainer.h"
40 #include "AliAODMCParticle.h"
41 #include "AliEMCALRecoUtils.h"
42 #include "AliAODEvent.h"
43 #include "AliESDEvent.h"
44 #include "AliAnalysisManager.h"
45 
47 
51 
52 // Actually registers the class with the base class
54 
55 const std::map <std::string, AliEMCALRecParam::AliEMCALClusterizerFlag> AliEmcalCorrectionClusterizer::fgkClusterizerTypeMap = {
56  {"kClusterizerv1", AliEMCALRecParam::kClusterizerv1 },
57  {"kClusterizerNxN", AliEMCALRecParam::kClusterizerNxN },
58  {"kClusterizerv2", AliEMCALRecParam::kClusterizerv2 },
59  {"kClusterizerv3", AliEMCALRecParam::kClusterizerv3 },
60  {"kClusterizerFW", AliEMCALRecParam::kClusterizerFW }
61 };
62 
67  AliEmcalCorrectionComponent("AliEmcalCorrectionClusterizer"),
68  fHistCPUTime(nullptr),
69  fHistRealTime(nullptr),
70  fTimer(nullptr),
71  fDigitsArr(0),
72  fClusterArr(0),
73  fRecParam(new AliEMCALRecParam),
74  fClusterizer(0),
75  fUnfolder(0),
76  fJustUnfold(kFALSE),
77  fUnfoldCellMinE(0),
78  fUnfoldCellMinEFrac(0),
79  fGeomName(),
80  fGeomMatrixSet(kFALSE),
81  fLoadGeomMatrices(kFALSE),
82  fOCDBpath(),
83  fCalibData(0),
84  fPedestalData(0),
85  fLoadCalib(kFALSE),
86  fLoadPed(kFALSE),
87  fSubBackground(kFALSE),
88  fNPhi(4),
89  fNEta(4),
90  fShiftPhi(2),
91  fShiftEta(2),
92  fTRUShift(0),
93  fTestPatternInput(kFALSE),
94  fSetCellMCLabelFromCluster(0),
95  fSetCellMCLabelFromEdepFrac(0),
96  fRemapMCLabelForAODs(0),
97  fRecalDistToBadChannels(kFALSE),
98  fRecalShowerShape(kFALSE),
99  fCaloClusters(0),
100  fEsd(0),
101  fAod(0)
102 {
103  for(Int_t i = 0; i < AliEMCALGeoParams::fgkEMCALModules; i++) fGeomMatrix[i] = 0 ;
104  for(Int_t j = 0; j < fgkTotalCellNumber; j++)
105  { fOrgClusterCellId[j] =-1; fCellLabels[j] =-1 ; }
106 }
107 
112 {
113  delete fClusterizer;
114  delete fUnfolder;
115  delete fRecParam;
116 }
117 
122 {
123  // Initialization
125 
126  std::string clusterizerTypeStr = "";
127  GetProperty("clusterizer", clusterizerTypeStr);
128  UInt_t clusterizerType = fgkClusterizerTypeMap.at(clusterizerTypeStr);
129 
130  Bool_t unfold = kFALSE;
131  GetProperty("unfold", unfold);
132  Bool_t unfoldRejectBelowThreshold = kFALSE;
133  GetProperty("unfoldRejectBelowThreshold", unfoldRejectBelowThreshold);
134  GetProperty("unfoldMinCellE" , fUnfoldCellMinE);
135  GetProperty("unfoldMinCellEFrac", fUnfoldCellMinEFrac);
136 
137  Double_t cellE = 0.05;
138  GetProperty("cellE", cellE);
139  Double_t seedE = 0.1;
140  GetProperty("seedE", seedE);
141  Float_t timeMin = -1; //minimum time of physical signal in a cell/digit (s) (in run macro, -50e-6)
142  GetProperty("cellTimeMin", timeMin);
143  Float_t timeMax = +1; //maximum time of physical signal in a cell/digit (s) (in run macro, 50e-6)
144  GetProperty("cellTimeMax", timeMax);
145  Float_t timeCut = 1; //maximum time difference between the digits inside EMC cluster (s) (in run macro, 1e-6)
146  GetProperty("clusterTimeLength", timeCut);
147  Float_t w0 = 4.5;
148  GetProperty("w0", w0);
149  GetProperty("recalDistToBadChannels", fRecalDistToBadChannels);
150  GetProperty("recalShowerShape", fRecalShowerShape);
151  GetProperty("remapMcAod", fRemapMCLabelForAODs);
152  Bool_t enableFracEMCRecalc = kFALSE;
153  GetProperty("enableFracEMCRecalc", enableFracEMCRecalc);
154  GetProperty("setCellMCLabelFromCluster", fSetCellMCLabelFromCluster);
155  Float_t diffEAggregation = 0.;
156  GetProperty("diffEAggregation", diffEAggregation);
157  GetProperty("useTestPatternForInput", fTestPatternInput);
158 
159  Int_t removeNMCGenerators = 0;
160  GetProperty("removeNMCGenerators", removeNMCGenerators);
161  Bool_t enableMCGenRemovTrack = 1;
162  GetProperty("enableMCGenRemovTrack", enableMCGenRemovTrack);
163  std::string removeMcGen1 = "";
164  GetProperty("removeMCGen1", removeMcGen1);
165  TString removeMCGen1 = removeMcGen1.c_str();
166  std::string removeMcGen2 = "";
167  GetProperty("removeMCGen2", removeMcGen2);
168  TString removeMCGen2 = removeMcGen2.c_str();
169 
170  fRecParam->SetClusterizerFlag(clusterizerType);
171  fRecParam->SetUnfold(unfold);
172  fRecParam->SetRejectBelowThreshold(unfoldRejectBelowThreshold);
173  fRecParam->SetMinECut(cellE);
174  fRecParam->SetClusteringThreshold(seedE);
175  fRecParam->SetW0(w0);
176  fRecParam->SetTimeMin(timeMin);
177  fRecParam->SetTimeMax(timeMax);
178  fRecParam->SetTimeCut(timeCut);
179  fRecParam->SetLocMaxCut(diffEAggregation); // Set minimum energy difference to start new cluster
180 
181  if (clusterizerType == AliEMCALRecParam::kClusterizerNxN)
182  fRecParam->SetNxM(1,1); // -> (1,1) means 3x3!
183 
184  // init reco utils
185  if (!fRecoUtils)
187 
188  if(enableFracEMCRecalc){
191 
192  printf("Enable frac MC recalc, remove generators %d \n",removeNMCGenerators);
193  if(removeNMCGenerators > 0)
194  {
195  printf("\t gen1 <%s>, gen2 <%s>, remove tracks %d\n",removeMCGen1.Data(),removeMCGen2.Data(),enableMCGenRemovTrack);
196  fRecoUtils->SetNumberOfMCGeneratorsToAccept(removeNMCGenerators) ;
197  fRecoUtils->SetNameOfMCGeneratorsToAccept(0,removeMCGen1);
198  fRecoUtils->SetNameOfMCGeneratorsToAccept(1,removeMCGen2);
199 
200  if(!enableMCGenRemovTrack) fRecoUtils->SwitchOffMCGeneratorToAcceptForTrackMatching();
201  }
202  }
203 
204  // Only support one cluster container for the clusterizer!
205  if (fClusterCollArray.GetEntries() > 1) {
206  AliFatal("Passed more than one cluster container to the clusterizer, but the clusterizer only supports one cluster container!");
207  }
208 
209  // Load 1D bad channel map
210  GetProperty("load1DBadChMap", fLoad1DBadChMap);
212 
213  return kTRUE;
214 }
215 
220 {
222 
223  if (fCreateHisto){
224  fHistCPUTime = new TH1F("hCPUTime","hCPUTime;CPU Time (ms)", 2000, 0, 1000);
225  fOutput->Add(fHistCPUTime);
226 
227  fHistRealTime = new TH1F("hRealTime","hRealTime;Real Time (ms)", 2000, 0, 1000);
228  fOutput->Add(fHistRealTime);
229 
230  fTimer = new TStopwatch();
231  }
232 }
233 
238 {
239  // Time the event loop if histogram creation is enabled
240  if (fCreateHisto) {
241  fTimer->Start(kTRUE);
242  }
243 
245 
246  fEsd = dynamic_cast<AliESDEvent*>(fEventManager.InputEvent());
247  fAod = dynamic_cast<AliAODEvent*>(fEventManager.InputEvent());
248 
249  // Only support one cluster container in the clusterizer!
251  if (!clusCont) {
252  AliFatal("Could not retrieve cluster container!");
253  }
254 
255  fCaloClusters = clusCont->GetArray();
256 
257  // If cells are empty, clear clusters and return
258  if (fCaloCells->GetNumberOfCells()<=0)
259  {
260  AliWarning(Form("Number of EMCAL cells = %d, returning", fCaloCells->GetNumberOfCells()));
262  return kFALSE;
263  }
264 
265  UInt_t offtrigger = 0;
266  if (fEsd) {
267  UInt_t mask1 = fEsd->GetESDRun()->GetDetectorsInDAQ();
268  UInt_t mask2 = fEsd->GetESDRun()->GetDetectorsInReco();
269  Bool_t desc1 = (mask1 >> 18) & 0x1;
270  Bool_t desc2 = (mask2 >> 18) & 0x1;
271  if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
272  AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
273  mask1, fEsd->GetESDRun()->GetDetectorsInReco(),
274  mask2, fEsd->GetESDRun()->GetDetectorsInDAQ()));
275  return kFALSE;
276  }
277  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
278  offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
279  } else {
280  offtrigger = ((AliVAODHeader*)fAod->GetHeader())->GetOfflineTrigger();
281  }
282 
283  if (!fMCEvent) {
284  if (offtrigger & AliVEvent::kFastOnly) {
285  AliError(Form("EMCAL not in fast only partition"));
286  return kFALSE;
287  }
288  }
289 
290  Init();
291 
292  if (fJustUnfold) {
293  AliWarning("Unfolding not implemented");
294  return kTRUE;
295  }
296 
297  FillDigitsArray();
298 
299  Clusterize();
300 
301  UpdateClusters();
302 
304 
305  if (fCreateHisto) {
306  fTimer->Stop();
307  // Ensure that it is stored in ms
308  fHistCPUTime->Fill(fTimer->CpuTime() * 1000.);
309  fHistRealTime->Fill(fTimer->RealTime() * 1000.);
310  }
311 
312  return kTRUE;
313 }
314 
319 {
320  if (fSubBackground) {
321  fClusterizer->SetInputCalibrated(kTRUE);
322  fClusterizer->SetCalibrationParameters(0);
323  }
324 
325  fClusterizer->Digits2Clusters("");
326 
327  if (fSubBackground) {
328  if (fCalibData) {
329  fClusterizer->SetInputCalibrated(kFALSE);
330  fClusterizer->SetCalibrationParameters(fCalibData);
331  }
332  }
333 }
334 
339 {
340  fDigitsArr->Clear("C");
341  if (fTestPatternInput) {
342  // Fill digits from a pattern
343  Int_t maxd = fGeom->GetNCells() / 4;
344  for (Int_t idigit = 0; idigit < maxd; idigit++){
345  if (idigit % 24 == 12) idigit += 12;
346  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
347  digit->SetId(idigit * 4);
348  digit->SetTime(600);
349  digit->SetTimeR(600);
350  digit->SetIndexInList(idigit);
351  digit->SetType(AliEMCALDigit::kHG);
352  digit->SetAmplitude(0.1);
353  }
354  }
355  else
356  {
357  // In case of MC productions done before aliroot tag v5-02-Rev09
358  // passing the cluster label to all the cells belonging to this cluster
359  // very rough
360  // Copied and simplified from AliEMCALTenderSupply
362  {
363  for (Int_t i = 0; i < fgkTotalCellNumber; i++)
364  {
365  fCellLabels [i] =-1 ;
366  fOrgClusterCellId[i] =-1 ;
367  }
368 
369  Int_t nClusters = fEventManager.InputEvent()->GetNumberOfCaloClusters();
370  for (Int_t i = 0; i < nClusters; i++)
371  {
372  AliVCluster *clus = fEventManager.InputEvent()->GetCaloCluster(i);
373 
374  if (!clus) continue;
375 
376  if (!clus->IsEMCAL()) continue ;
377 
378  Int_t label = clus->GetLabel();
379  UShort_t * index = clus->GetCellsAbsId() ;
380 
381  for(Int_t icell=0; icell < clus->GetNCells(); icell++)
382  {
384  fCellLabels[index[icell]] = label;
385 
386  fOrgClusterCellId[index[icell]] = i ; // index of the original cluster
387  } // cell in cluster loop
388  } // cluster loop
389  }
390 
391  Double_t avgE = 0; // for background subtraction
392  const Int_t ncells = fCaloCells->GetNumberOfCells();
393  for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell)
394  {
395  Float_t cellAmplitude=0;
396  Double_t cellTime=0, amp = 0, cellEFrac = 0;
397  Short_t cellNumber=0;
398  Int_t cellMCLabel=-1;
399  if (fCaloCells->GetCell(icell, cellNumber, amp, cellTime, cellMCLabel, cellEFrac) != kTRUE)
400  break;
401 
402  cellAmplitude = amp; // compilation problem
403 
404  //if (fSetCellMCLabelFromCluster) cellMCLabel = fCellLabels[cellNumber];
406  {
407  if (fSetCellMCLabelFromCluster) cellMCLabel = fCellLabels[cellNumber];
408  else if (fRemapMCLabelForAODs ) RemapMCLabelForAODs(cellMCLabel);
409  }
410 
411  if (cellMCLabel > 0 && cellEFrac < 1e-6)
412  cellEFrac = 1;
413 
414  if (cellAmplitude < 1e-6 || cellNumber < 0)
415  continue;
416 
417  // New way to set the cell MC labels,
418  // valid only for MC productions with aliroot > v5-07-21
419  //
420  TArrayI labeArr(0);
421  TArrayF eDepArr(0);
422  Int_t nLabels = 0;
423 
424  if(fSetCellMCLabelFromEdepFrac && fOrgClusterCellId[cellNumber] >= 0) // index can be negative if noisy cell that did not form cluster
425  {
426  cellMCLabel = -1;
427 
428  fCellLabels[cellNumber] = idigit;
429 
430  Int_t iclus = fOrgClusterCellId[cellNumber];
431 
432  if(iclus < 0)
433  {
434  AliInfo("Negative original cluster index, skip \n");
435  continue;
436  }
437 
438  AliVCluster* clus = fEventManager.InputEvent()->GetCaloCluster(iclus);
439 
440  for(Int_t icluscell=0; icluscell < clus->GetNCells(); icluscell++ )
441  {
442  if(cellNumber != clus->GetCellAbsId(icluscell)) continue ;
443 
444  // Select the MC label contributing, only if enough energy deposition
445  fRecoUtils->RecalculateCellLabelsRemoveAddedGenerator(cellNumber, clus, fMCEvent, cellAmplitude, labeArr, eDepArr);
446  nLabels = labeArr.GetSize();
447  }
448  }
449 
450  AliEMCALDigit *digit = new((*fDigitsArr)[idigit]) AliEMCALDigit(cellMCLabel, cellMCLabel, cellNumber,
451  cellAmplitude, (Float_t)cellTime,
452  AliEMCALDigit::kHG,idigit, 0, 0, cellEFrac*cellAmplitude);
453 
454  if(nLabels > 0)
455  {
456  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
457  }
458 
459  if (fSubBackground)
460  {
461  Float_t energy = cellAmplitude;
462  Float_t time = cellTime;
463  fClusterizer->Calibrate(energy,time,cellNumber);
464  digit->SetAmplitude(energy);
465  avgE += energy;
466  }
467 
468  idigit++;
469  }
470 
471  if (fSubBackground) {
472  avgE /= fGeom->GetNumberOfSuperModules()*48*24;
473  Int_t ndigis = fDigitsArr->GetEntries();
474  for (Int_t i = 0; i < ndigis; ++i) {
475  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
476  Double_t energy = digit->GetAmplitude() - avgE;
477  if (energy<=0.001) {
478  digit->SetAmplitude(0);
479  } else {
480  digit->SetAmplitude(energy);
481  }
482  }
483  }
484  }
485 }
486 
492 {
493  const Int_t Ncls = fClusterArr->GetEntries();
494  AliDebug(1, Form("total no of clusters %d", Ncls));
495 
496  for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i)
497  {
498  AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
499 
500  Int_t ncellsTrue = 0;
501  const Int_t ncells = recpoint->GetMultiplicity();
502  UShort_t absIds[ncells];
503  Double32_t ratios[ncells];
504  Int_t *dlist = recpoint->GetDigitsList();
505  Float_t *elist = recpoint->GetEnergiesList();
506  Double_t mcEnergy = 0;
507 
508  Float_t clusterE = 0;
509  for (Int_t c = 0; c < ncells; ++c)
510  {
511  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
512  absIds[ncellsTrue] = digit->GetId();
513  ratios[ncellsTrue] = elist[c]/digit->GetAmplitude();
514 
515  if ( !fRecParam->GetUnfold() && (ratios[ncellsTrue] > 1 || ratios[ncellsTrue] < 1) )
516  AliWarning(Form("recpoint cell E %2.3f but digit E %2.3f and no unfolding",
517  recpoint->GetEnergiesList()[c], digit->GetAmplitude()));
518 
519  // In case of unfolding, remove digits with unfolded energy too low
520  if ( fRecParam->GetUnfold() )
521  {
522  if ( recpoint->GetEnergiesList()[c] < fUnfoldCellMinE ||
523  ratios[ncellsTrue] < fUnfoldCellMinEFrac )
524  {
525 
526  AliDebug(2,Form("Too small energy in cell of cluster: cluster cell %f, digit %f",
527  recpoint->GetEnergiesList()[c],digit->GetAmplitude()));
528  continue;
529  } // if cuts
530  }// Select cells
531 
532  // Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
533  clusterE += recpoint->GetEnergiesList()[c];
534 
535  if (digit->GetIparent(1) > 0)
536  mcEnergy += digit->GetDEParent(1)/recpoint->GetEnergy();
537 
538  ++ncellsTrue;
539  } // cluster cell loop
540 
541  if (ncellsTrue < 1)
542  {
543  AliWarning("Skipping cluster with no cells");
544  continue;
545  }
546 
547  // calculate new cluster position
548  TVector3 gpos;
549  recpoint->GetGlobalPosition(gpos);
550  Float_t g[3];
551  gpos.GetXYZ(g);
552 
553  AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
554 
555  AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
556  c->SetType(AliVCluster::kEMCALClusterv1);
557  c->SetE(clusterE);
558  c->SetPosition(g);
559  c->SetNCells(ncellsTrue);
560  c->SetCellsAbsId(absIds);
561  c->SetCellsAmplitudeFraction(ratios);
562  c->SetID(nout-1);
563  c->SetDispersion(recpoint->GetDispersion());
564  c->SetEmcCpvDistance(-1);
565  c->SetChi2(-1);
566  c->SetTOF(recpoint->GetTime()) ; //time-of-flight
567  c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
568  c->SetMCEnergyFraction(mcEnergy);
569 
570  if ( ncells == ncellsTrue )
571  {
572  Float_t elipAxis[2];
573  recpoint->GetElipsAxis(elipAxis);
574  c->SetM02(elipAxis[0]*elipAxis[0]);
575  c->SetM20(elipAxis[1]*elipAxis[1]);
576  }
577  else
578  {
579  // In case some cells rejected, in unfolding case, recalculate
580  // shower shape parameters and position
581  AliDebug(2,Form("Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape",ncells,ncellsTrue));
582 
584  //fRecoUtils->RecalculateClusterPID(c);
586  }
587 
588  //
589  // MC labels
590  //
591  Int_t parentMult = 0;
592  Int_t *parentList = recpoint->GetParents(parentMult);
593  Float_t *parentListDE = recpoint->GetParentsDE(); // deposited energy
594 
595  c->SetLabel(parentList, parentMult);
596  c->SetClusterMCEdepFractionFromEdepArray(parentListDE);
597 
598  //
599  // Set the cell energy deposition fraction map:
600  //
601  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
602  {
603  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
604 
605  // Get the digit that originated this cell cluster
606  //AliVCaloCells* cells = InputEvent()->GetEMCALCells();
607 
608  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
609  {
610  Int_t idigit = fCellLabels[absIds[icell]];
611 
612  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
613 
614  // Find the 4 MC labels that contributed to the cluster and their
615  // deposited energy in the current digit
616 
617  mcEdepFracPerCell[icell] = 0; // init
618 
619  Int_t nparents = dig->GetNiparent();
620  if ( nparents > 0 )
621  {
622  Int_t digLabel =-1 ;
623  Float_t edep = 0 ;
624  Float_t edepTot = 0 ;
625  Float_t mcEDepFrac[4] = {0,0,0,0};
626 
627  // all parents in digit
628  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
629  {
630  digLabel = dig->GetIparent (jndex+1);
631  edep = dig->GetDEParent(jndex+1);
632  edepTot += edep;
633 
634  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
635  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
636  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
637  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
638  } // all prarents in digit
639 
640  // Divide energy deposit by total deposited energy
641  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
642  if(edepTot > 0.01)
643  {
644  mcEdepFracPerCell[icell] = c->PackMCEdepFraction(mcEDepFrac);
645  }
646  } // at least one parent label in digit
647  } // cell in cluster loop
648 
649  c->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
650 
651  delete [] mcEdepFracPerCell;
652 
653  } // at least one parent in cluster, do the cell primary packing
654  }
655 }
656 
661 {
662  // Before destroying the orignal list, assign to the rec points the MC labels
663  // of the original clusters, if requested
666 
668 
669  fCaloClusters->Compress();
670 
672 }
673 
678 {
679  Int_t nclusters = fCaloClusters->GetEntriesFast();
680  for (Int_t icluster=0; icluster < nclusters; ++icluster) {
681  AliVCluster *clust = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
682  if (!clust) {
683  continue;
684  }
685  if (!clust->IsEMCAL()) {
686  continue;
687  }
688 
689  // SHOWER SHAPE -----------------------------------------------
690  if (fRecalShowerShape)
692 
693  // DISTANCE TO BAD CHANNELS -----------------------------------
696 
697  }
698 
699  fCaloClusters->Compress();
700 }
701 
706 {
707  if (label < 0) return;
708 
709  TClonesArray * arr = dynamic_cast<TClonesArray*>(fAod->FindListObject("mcparticles")) ;
710  if (!arr) return ;
711 
712  if (label < arr->GetEntriesFast())
713  {
714  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
715  if (!particle) return ;
716 
717  if (label == particle->Label()) return ; // label already OK
718  }
719 
720  // loop on the particles list and check if there is one with the same label
721  for (Int_t ind = 0; ind < arr->GetEntriesFast(); ind++)
722  {
723  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
724  if (!particle) continue ;
725 
726  if (label == particle->Label())
727  {
728  label = ind;
729  return;
730  }
731  }
732 
733  label = -1;
734 }
735 
745 {
746  Int_t ncls = fClusterArr->GetEntriesFast();
747  for(Int_t irp=0; irp < ncls; ++irp)
748  {
749  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
750  clArray.Reset();
751  Int_t nClu = 0;
752  Int_t nLabTotOrg = 0;
753  Float_t emax = -1;
754  Int_t idMax = -1;
755 
756  AliEMCALRecPoint *clus = static_cast<AliEMCALRecPoint*>(fClusterArr->At(irp));
757 
758  //Find the clusters that originally had the cells
759  const Int_t ncells = clus->GetMultiplicity();
760  Int_t *digList = clus->GetDigitsList();
761 
762  for (Int_t iLoopCell = 0 ; iLoopCell < ncells ; iLoopCell++)
763  {
764  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(digList[iLoopCell]));
765  Int_t idCell = digit->GetId();
766 
767  if (idCell>=0)
768  {
769  Int_t idCluster = fOrgClusterCellId[idCell];
770  Bool_t set = kTRUE;
771  for (Int_t icl =0; icl < nClu; icl++)
772  {
773  if (((Int_t)clArray.GetAt(icl))==-1) continue;
774  if (idCluster == ((Int_t)clArray.GetAt(icl))) set = kFALSE;
775  }
776  if (set && idCluster >= 0)
777  {
778  clArray.SetAt(idCluster,nClu++);
779  nLabTotOrg+=(fEventManager.InputEvent()->GetCaloCluster(idCluster))->GetNLabels();
780 
781  //Search highest E cluster
782  AliVCluster * clOrg = fEventManager.InputEvent()->GetCaloCluster(idCluster);
783  if (emax < clOrg->E())
784  {
785  emax = clOrg->E();
786  idMax = idCluster;
787  }
788  }
789  }
790  }// cell loop
791 
792  // Put the first in the list the cluster with highest energy
793  if (idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
794  {
795  Int_t maxIndex = -1;
796  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
797  for (Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++)
798  {
799  if (idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
800  }
801 
802  if (firstCluster >=0 && idMax >=0)
803  {
804  clArray.SetAt(idMax,0);
805  clArray.SetAt(firstCluster,maxIndex);
806  }
807  }
808 
809  // Get the labels list in the original clusters, assign all to the new cluster
810  TArrayI clMCArray(nLabTotOrg) ;
811  clMCArray.Reset();
812 
813  Int_t nLabTot = 0;
814  for (Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++)
815  {
816  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
817  AliVCluster * clOrg = fEventManager.InputEvent()->GetCaloCluster(idCluster);
818  Int_t nLab = clOrg->GetNLabels();
819 
820  for (Int_t iLab = 0 ; iLab < nLab ; iLab++)
821  {
822  Int_t lab = clOrg->GetLabelAt(iLab) ;
823  if (lab>=0)
824  {
825  Bool_t set = kTRUE;
826  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
827  {
828  if (lab == ((Int_t)clMCArray.GetAt(iLabTot))) set = kFALSE;
829  }
830  if (set) clMCArray.SetAt(lab,nLabTot++);
831  }
832  }
833  }// cluster loop
834 
835  // Set the final list of labels to rec point
836 
837  Int_t *labels = new Int_t[nLabTot];
838  for(Int_t il = 0; il < nLabTot; il++) labels[il] = clMCArray.GetArray()[il];
839  clus->SetParents(nLabTot,labels);
840 
841  } // rec point array
842 }
843 
848 {
849  // Select clusterization/unfolding algorithm and set all the needed parameters.
850 
851  // If distBC option enabled, init and fill bad channel map
854  else
856 
858 
859  if (fJustUnfold){
860  // init the unfolding afterburner
861  delete fUnfolder;
862  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
863  return;
864  }
865 
866  if (!fGeomMatrixSet) {
867  if (fLoadGeomMatrices) {
868  for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
869  if (fGeomMatrix[mod]){
870  if (AliAnalysisManager::GetAnalysisManager()->GetDebugLevel() > 2)
871  fGeomMatrix[mod]->Print();
872  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
873  }
874  }
875  } else { // get matrix from file (work around bug in aliroot)
876  for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
877  const TGeoHMatrix *gm = 0;
878  if (fEsd) {
879  gm = fEsd->GetEMCALMatrix(mod);
880  } else {
881  AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(fAod->GetHeader());
882  if(!aodheader) AliFatal("Not a standard AOD");
883  if (aodheader) {
884  gm = aodheader->GetEMCALMatrix(mod);
885  }
886  }
887  if (gm) {
888  if (AliAnalysisManager::GetAnalysisManager()->GetDebugLevel() > 2)
889  gm->Print();
890  fGeom->SetMisalMatrix(gm,mod);
891  }
892  }
893  }
894  fGeomMatrixSet=kTRUE;
895  }
896 
897  // setup digit array if needed
898  if (!fDigitsArr) {
899  fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
900  fDigitsArr->SetOwner(1);
901  }
902 
903  // then setup clusterizer
904  if (fClusterizer) {
905  // avoid to delete digits array
906  fClusterizer->SetDigitsArr(0);
907  delete fClusterizer;
908  }
909  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
910  fClusterizer = new AliEMCALClusterizerv1(fGeom);
911  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
912  AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
913  clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
914  clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
915  fClusterizer = clusterizer;
916  }
917  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
918  fClusterizer = new AliEMCALClusterizerv2(fGeom);
919  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv3)
920  fClusterizer = new AliEMCALClusterizerv3(fGeom);
921  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
922  AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
923  clusterizer->SetNphi(fNPhi);
924  clusterizer->SetNeta(fNEta);
925  clusterizer->SetShiftPhi(fShiftPhi);
926  clusterizer->SetShiftEta(fShiftEta);
927  clusterizer->SetTRUshift(fTRUShift);
928  fClusterizer = clusterizer;
929  }
930  else {
931  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
932  }
933  fClusterizer->InitParameters(fRecParam);
934 
935  if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
936  AliCDBManager *cdb = AliCDBManager::Instance();
937  if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
938  cdb->SetDefaultStorage(fOCDBpath);
939  if (fRun!=cdb->GetRun())
940  cdb->SetRun(fRun);
941  }
942  if (!fCalibData&&fLoadCalib&&fRun>0) {
943  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
944  if (entry)
945  fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
946  if (!fCalibData)
947  AliFatal("Calibration parameters not found in CDB!");
948  }
949  if (!fPedestalData&&fLoadPed&&fRun>0) {
950  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
951  if (entry)
952  fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
953  }
954  if (fCalibData) {
955  fClusterizer->SetInputCalibrated(kFALSE);
956  fClusterizer->SetCalibrationParameters(fCalibData);
957  } else {
958  fClusterizer->SetInputCalibrated(kTRUE);
959  }
960  fClusterizer->SetCaloCalibPedestal(fPedestalData);
961  fClusterizer->SetJustClusters(kTRUE);
962  fClusterizer->SetDigitsArr(fDigitsArr);
963  fClusterizer->SetOutput(0);
964  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
965 
966  //fClusterizer->Print("");
967 }
968 
974 {
976 
977  if (runChanged && fRecalDistToBadChannels) {
978  // init bad channels
979  Int_t fInitBC = InitBadChannels();
980  if (fInitBC==0) {
981  AliError("InitBadChannels returned false, returning");
982  }
983  if (fInitBC==1) {
984  AliWarning("InitBadChannels OK");
985  }
986  if (fInitBC>1) {
987  AliWarning(Form("No external hot channel set: %d - %s", fEventManager.InputEvent()->GetRunNumber(), fFilepass.Data()));
988  }
989  }
990  return runChanged;
991 }
992 
997 {
998  const Int_t nents = fCaloClusters->GetEntries();
999  for (Int_t i=0;i<nents;++i) {
1000  AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
1001  if (!c) {
1002  continue;
1003  }
1004  if (c->IsEMCAL()) {
1005  delete fCaloClusters->RemoveAt(i);
1006  }
1007  }
1008 }
Bool_t fTestPatternInput
Use test pattern as input instead of cells.
void RecalculateCellLabelsRemoveAddedGenerator(Int_t absID, AliVCluster *clus, AliMCEvent *mc, Float_t &amp, TArrayI &labeArr, TArrayF &eDepArr) const
TStopwatch * fTimer
! Timer for the Run() function (event loop)
Bool_t fRemapMCLabelForAODs
Remap AOD cells MC label.
Int_t fShiftPhi
shift in phi (for FixedWindowsClusterizer)
void Print(std::ostream &o, const char *name, Double_t dT, Double_t dVM, Double_t alldT, Double_t alldVM)
Definition: PlotSysInfo.C:121
AliEMCALGeometry * fGeom
! Geometry object
static const Int_t fgkTotalCellNumber
Maximum number of cells in EMCAL/DCAL: (48*24)*(10+4/3.+6*2/3.)
double Double_t
Definition: External.C:58
TClonesArray * fCaloClusters
!calo clusters array
AliEMCALRecParam * fRecParam
reconstruction parameters container
Bool_t fSetCellMCLabelFromEdepFrac
For MC generated with aliroot > v5-07-21, check the EDep information.
Int_t fCellLabels[fgkTotalCellNumber]
Array with MC label/map.
void SwitchOffDistToBadChannelRecalculation()
AliEMCALAfterBurnerUF * fUnfolder
!unfolding procedure
Bool_t fLoadPed
access ped object from OCDB (def=off)
Bool_t fTRUShift
shifting inside a TRU (true) or through the whole calorimeter (false) (for FixedWindowsClusterizer) ...
energy
Definition: HFPtSpectrum.C:45
void SetNumberOfMCGeneratorsToAccept(Int_t nGen)
Int_t fShiftEta
shift in eta (for FixedWindowsClusterizer)
TObjArray * fClusterArr
!recpoints array
TCanvas * c
Definition: TestFitELoss.C:172
Int_t fNPhi
nPhi (for FixedWindowsClusterizer)
AliVCaloCells * fCaloCells
! Pointer to CaloCells
Some utilities for cluster and cell treatment.
void SwitchOffMCGeneratorToAcceptForTrackMatching()
Bool_t fLoadCalib
access calib object from OCDB (def=off)
AliEMCALRecoUtils * fRecoUtils
Pointer to RecoUtils.
void RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
TObjArray fClusterCollArray
Cluster collection array.
AliCaloCalibPedestal * fPedestalData
EMCAL pedestal.
TGeoHMatrix * fGeomMatrix[AliEMCALGeoParams::fgkEMCALModules]
geometry matrices with alignments
AliEMCALClusterizer * fClusterizer
!clusterizer
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
void RecalculateClusterShowerShapeParameters(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
float Float_t
Definition: External.C:68
Base class for correction components in the EMCal correction framework.
EMCal clusterizer component in the EMCal correction framework.
Int_t fSetCellMCLabelFromCluster
Use cluster MC label as cell label:
Bool_t fRecalShowerShape
switch for recalculation of the shower shape
TClonesArray * GetArray() const
Bool_t fSubBackground
subtract background if true (def=off)
Int_t fNEta
nEta (for FixedWinoswsClusterizer)
TList * fOutput
! List of output histograms
short Short_t
Definition: External.C:23
void SwitchOnDistToBadChannelRecalculation()
Bool_t fGeomMatrixSet
set geometry matrices only once, for the first event.
Bool_t fCreateHisto
Flag to make some basic histograms.
TString fOCDBpath
path with OCDB location
Float_t fUnfoldCellMinE
min energy cell threshold, after unfolding
static RegisterCorrectionComponent< AliEmcalCorrectionClusterizer > reg
Bool_t fRecalDistToBadChannels
recalculate distance to bad channel
static const std::map< std::string, AliEMCALRecParam::AliEMCALClusterizerFlag > fgkClusterizerTypeMap
Relates string to the clusterizer type enumeration for YAML configuration.
void RecalculateClusterPosition(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
TH1F * fHistCPUTime
! CPU time for the Run() function (event loop)
AliEMCALCalibData * fCalibData
EMCAL calib data.
Int_t fOrgClusterCellId[fgkTotalCellNumber]
Array ID of cluster to which the cell belongs in unmodified clusters.
AliEmcalCorrectionEventManager fEventManager
Minimal task which inherits from AliAnalysisTaskSE and manages access to the event.
void SetUse1DBadChannelMap(Bool_t use)
Float_t fUnfoldCellMinEFrac
min fraction of cell energy after unfolding
Bool_t fJustUnfold
just unfold, do not recluster
unsigned short UShort_t
Definition: External.C:28
Bool_t fLoadGeomMatrices
matrices from configuration, not geometry.root nor ESDs/AODs
bool Bool_t
Definition: External.C:53
Bool_t fLoad1DBadChMap
Flag to load 1D bad channel map.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
TH1F * fHistRealTime
! Real time for the Run() function (event loop)
Container structure for EMCAL clusters.
TClonesArray * fDigitsArr
!digits array
TString fFilepass
Input data pass number.
bool GetProperty(std::string propertyName, T &property, bool requiredProperty=true, std::string correctionName="")
Retrieve property.
void SetNameOfMCGeneratorsToAccept(Int_t ig, TString name)