AliPhysics  64f4410 (64f4410)
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  fGeomName(),
78  fGeomMatrixSet(kFALSE),
79  fLoadGeomMatrices(kFALSE),
80  fOCDBpath(),
81  fCalibData(0),
82  fPedestalData(0),
83  fLoadCalib(kFALSE),
84  fLoadPed(kFALSE),
85  fSubBackground(kFALSE),
86  fNPhi(4),
87  fNEta(4),
88  fShiftPhi(2),
89  fShiftEta(2),
90  fTRUShift(0),
91  fTestPatternInput(kFALSE),
92  fSetCellMCLabelFromCluster(0),
93  fSetCellMCLabelFromEdepFrac(0),
94  fRemapMCLabelForAODs(0),
95  fRecalDistToBadChannels(kFALSE),
96  fRecalShowerShape(kFALSE),
97  fCaloClusters(0),
98  fEsd(0),
99  fAod(0)
100 {
101  for(Int_t i = 0; i < AliEMCALGeoParams::fgkEMCALModules; i++) fGeomMatrix[i] = 0 ;
102  for(Int_t j = 0; j < fgkTotalCellNumber; j++)
103  { fOrgClusterCellId[j] =-1; fCellLabels[j] =-1 ; }
104 }
105 
110 {
111  delete fClusterizer;
112  delete fUnfolder;
113  delete fRecParam;
114 }
115 
120 {
121  // Initialization
123 
124  std::string clusterizerTypeStr = "";
125  GetProperty("clusterizer", clusterizerTypeStr);
126  UInt_t clusterizerType = fgkClusterizerTypeMap.at(clusterizerTypeStr);
127  Double_t cellE = 0.05;
128  GetProperty("cellE", cellE);
129  Double_t seedE = 0.1;
130  GetProperty("seedE", seedE);
131  Float_t timeMin = -1; //minimum time of physical signal in a cell/digit (s) (in run macro, -50e-6)
132  GetProperty("cellTimeMin", timeMin);
133  Float_t timeMax = +1; //maximum time of physical signal in a cell/digit (s) (in run macro, 50e-6)
134  GetProperty("cellTimeMax", timeMax);
135  Float_t timeCut = 1; //maximum time difference between the digits inside EMC cluster (s) (in run macro, 1e-6)
136  GetProperty("clusterTimeLength", timeCut);
137  Float_t w0 = 4.5;
138  GetProperty("w0", w0);
139  GetProperty("recalDistToBadChannels", fRecalDistToBadChannels);
140  GetProperty("recalShowerShape", fRecalShowerShape);
141  GetProperty("remapMcAod", fRemapMCLabelForAODs);
142  Bool_t enableFracEMCRecalc = kFALSE;
143  GetProperty("enableFracEMCRecalc", enableFracEMCRecalc);
144  GetProperty("setCellMCLabelFromCluster", fSetCellMCLabelFromCluster);
145  Float_t diffEAggregation = 0.;
146  GetProperty("diffEAggregation", diffEAggregation);
147  GetProperty("useTestPatternForInput", fTestPatternInput);
148 
149  Int_t removeNMCGenerators = 0;
150  GetProperty("removeNMCGenerators", removeNMCGenerators);
151  Bool_t enableMCGenRemovTrack = 1;
152  GetProperty("enableMCGenRemovTrack", enableMCGenRemovTrack);
153  std::string removeMcGen1 = "";
154  GetProperty("removeMCGen1", removeMcGen1);
155  TString removeMCGen1 = removeMcGen1.c_str();
156  std::string removeMcGen2 = "";
157  GetProperty("removeMCGen2", removeMcGen2);
158  TString removeMCGen2 = removeMcGen2.c_str();
159 
160  fRecParam->SetClusterizerFlag(clusterizerType);
161  fRecParam->SetMinECut(cellE);
162  fRecParam->SetClusteringThreshold(seedE);
163  fRecParam->SetW0(w0);
164  fRecParam->SetTimeMin(timeMin);
165  fRecParam->SetTimeMax(timeMax);
166  fRecParam->SetTimeCut(timeCut);
167  fRecParam->SetLocMaxCut(diffEAggregation); // Set minimum energy difference to start new cluster
168 
169  if (clusterizerType == AliEMCALRecParam::kClusterizerNxN)
170  fRecParam->SetNxM(1,1); // -> (1,1) means 3x3!
171 
172  // init reco utils
173  if (!fRecoUtils)
175 
176  if(enableFracEMCRecalc){
179 
180  printf("Enable frac MC recalc, remove generators %d \n",removeNMCGenerators);
181  if(removeNMCGenerators > 0)
182  {
183  printf("\t gen1 <%s>, gen2 <%s>, remove tracks %d\n",removeMCGen1.Data(),removeMCGen2.Data(),enableMCGenRemovTrack);
184  fRecoUtils->SetNumberOfMCGeneratorsToAccept(removeNMCGenerators) ;
185  fRecoUtils->SetNameOfMCGeneratorsToAccept(0,removeMCGen1);
186  fRecoUtils->SetNameOfMCGeneratorsToAccept(1,removeMCGen2);
187 
188  if(!enableMCGenRemovTrack) fRecoUtils->SwitchOffMCGeneratorToAcceptForTrackMatching();
189  }
190  }
191 
192  // Only support one cluster container for the clusterizer!
193  if (fClusterCollArray.GetEntries() > 1) {
194  AliFatal("Passed more than one cluster container to the clusterizer, but the clusterizer only supports one cluster container!");
195  }
196 
197  return kTRUE;
198 }
199 
204 {
206 
207  if (fCreateHisto){
208  fHistCPUTime = new TH1F("hCPUTime","hCPUTime;CPU Time (ms)", 2000, 0, 1000);
209  fOutput->Add(fHistCPUTime);
210 
211  fHistRealTime = new TH1F("hRealTime","hRealTime;Real Time (ms)", 2000, 0, 1000);
212  fOutput->Add(fHistRealTime);
213 
214  fTimer = new TStopwatch();
215  }
216 }
217 
222 {
223  // Time the event loop if histogram creation is enabled
224  if (fCreateHisto) {
225  fTimer->Start(kTRUE);
226  }
227 
229 
230  fEsd = dynamic_cast<AliESDEvent*>(fEventManager.InputEvent());
231  fAod = dynamic_cast<AliAODEvent*>(fEventManager.InputEvent());
232 
233  // Only support one cluster container in the clusterizer!
235  if (!clusCont) {
236  AliFatal("Could not retrieve cluster container!");
237  }
238 
239  fCaloClusters = clusCont->GetArray();
240 
241  // If cells are empty, clear clusters and return
242  if (fCaloCells->GetNumberOfCells()<=0)
243  {
244  AliWarning(Form("Number of EMCAL cells = %d, returning", fCaloCells->GetNumberOfCells()));
246  return kFALSE;
247  }
248 
249  UInt_t offtrigger = 0;
250  if (fEsd) {
251  UInt_t mask1 = fEsd->GetESDRun()->GetDetectorsInDAQ();
252  UInt_t mask2 = fEsd->GetESDRun()->GetDetectorsInReco();
253  Bool_t desc1 = (mask1 >> 18) & 0x1;
254  Bool_t desc2 = (mask2 >> 18) & 0x1;
255  if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
256  AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
257  mask1, fEsd->GetESDRun()->GetDetectorsInReco(),
258  mask2, fEsd->GetESDRun()->GetDetectorsInDAQ()));
259  return kFALSE;
260  }
261  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
262  offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
263  } else {
264  offtrigger = ((AliVAODHeader*)fAod->GetHeader())->GetOfflineTrigger();
265  }
266 
267  if (!fMCEvent) {
268  if (offtrigger & AliVEvent::kFastOnly) {
269  AliError(Form("EMCAL not in fast only partition"));
270  return kFALSE;
271  }
272  }
273 
274  Init();
275 
276  if (fJustUnfold) {
277  AliWarning("Unfolding not implemented");
278  return kTRUE;
279  }
280 
281  FillDigitsArray();
282 
283  Clusterize();
284 
285  UpdateClusters();
286 
288 
289  if (fCreateHisto) {
290  fTimer->Stop();
291  // Ensure that it is stored in ms
292  fHistCPUTime->Fill(fTimer->CpuTime() * 1000.);
293  fHistRealTime->Fill(fTimer->RealTime() * 1000.);
294  }
295 
296  return kTRUE;
297 }
298 
303 {
304  if (fSubBackground) {
305  fClusterizer->SetInputCalibrated(kTRUE);
306  fClusterizer->SetCalibrationParameters(0);
307  }
308 
309  fClusterizer->Digits2Clusters("");
310 
311  if (fSubBackground) {
312  if (fCalibData) {
313  fClusterizer->SetInputCalibrated(kFALSE);
314  fClusterizer->SetCalibrationParameters(fCalibData);
315  }
316  }
317 }
318 
323 {
324  fDigitsArr->Clear("C");
325  if (fTestPatternInput) {
326  // Fill digits from a pattern
327  Int_t maxd = fGeom->GetNCells() / 4;
328  for (Int_t idigit = 0; idigit < maxd; idigit++){
329  if (idigit % 24 == 12) idigit += 12;
330  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
331  digit->SetId(idigit * 4);
332  digit->SetTime(600);
333  digit->SetTimeR(600);
334  digit->SetIndexInList(idigit);
335  digit->SetType(AliEMCALDigit::kHG);
336  digit->SetAmplitude(0.1);
337  }
338  }
339  else
340  {
341  // In case of MC productions done before aliroot tag v5-02-Rev09
342  // passing the cluster label to all the cells belonging to this cluster
343  // very rough
344  // Copied and simplified from AliEMCALTenderSupply
346  {
347  for (Int_t i = 0; i < fgkTotalCellNumber; i++)
348  {
349  fCellLabels [i] =-1 ;
350  fOrgClusterCellId[i] =-1 ;
351  }
352 
353  Int_t nClusters = fEventManager.InputEvent()->GetNumberOfCaloClusters();
354  for (Int_t i = 0; i < nClusters; i++)
355  {
356  AliVCluster *clus = fEventManager.InputEvent()->GetCaloCluster(i);
357 
358  if (!clus) continue;
359 
360  if (!clus->IsEMCAL()) continue ;
361 
362  Int_t label = clus->GetLabel();
363  UShort_t * index = clus->GetCellsAbsId() ;
364 
365  for(Int_t icell=0; icell < clus->GetNCells(); icell++)
366  {
368  fCellLabels[index[icell]] = label;
369 
370  fOrgClusterCellId[index[icell]] = i ; // index of the original cluster
371  } // cell in cluster loop
372  } // cluster loop
373  }
374 
375  Double_t avgE = 0; // for background subtraction
376  const Int_t ncells = fCaloCells->GetNumberOfCells();
377  for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell)
378  {
379  Float_t cellAmplitude=0;
380  Double_t cellTime=0, amp = 0, cellEFrac = 0;
381  Short_t cellNumber=0;
382  Int_t cellMCLabel=-1;
383  if (fCaloCells->GetCell(icell, cellNumber, amp, cellTime, cellMCLabel, cellEFrac) != kTRUE)
384  break;
385 
386  cellAmplitude = amp; // compilation problem
387 
388  //if (fSetCellMCLabelFromCluster) cellMCLabel = fCellLabels[cellNumber];
390  {
391  if (fSetCellMCLabelFromCluster) cellMCLabel = fCellLabels[cellNumber];
392  else if (fRemapMCLabelForAODs ) RemapMCLabelForAODs(cellMCLabel);
393  }
394 
395  if (cellMCLabel > 0 && cellEFrac < 1e-6)
396  cellEFrac = 1;
397 
398  if (cellAmplitude < 1e-6 || cellNumber < 0)
399  continue;
400 
401  // New way to set the cell MC labels,
402  // valid only for MC productions with aliroot > v5-07-21
403  //
404  TArrayI labeArr(0);
405  TArrayF eDepArr(0);
406  Int_t nLabels = 0;
407 
408  if(fSetCellMCLabelFromEdepFrac && fOrgClusterCellId[cellNumber] >= 0) // index can be negative if noisy cell that did not form cluster
409  {
410  cellMCLabel = -1;
411 
412  fCellLabels[cellNumber] = idigit;
413 
414  Int_t iclus = fOrgClusterCellId[cellNumber];
415 
416  if(iclus < 0)
417  {
418  AliInfo("Negative original cluster index, skip \n");
419  continue;
420  }
421 
422  AliVCluster* clus = fEventManager.InputEvent()->GetCaloCluster(iclus);
423 
424  for(Int_t icluscell=0; icluscell < clus->GetNCells(); icluscell++ )
425  {
426  if(cellNumber != clus->GetCellAbsId(icluscell)) continue ;
427 
428  // Select the MC label contributing, only if enough energy deposition
429  fRecoUtils->RecalculateCellLabelsRemoveAddedGenerator(cellNumber, clus, fMCEvent, cellAmplitude, labeArr, eDepArr);
430  nLabels = labeArr.GetSize();
431  }
432  }
433 
434  AliEMCALDigit *digit = new((*fDigitsArr)[idigit]) AliEMCALDigit(cellMCLabel, cellMCLabel, cellNumber,
435  cellAmplitude, (Float_t)cellTime,
436  AliEMCALDigit::kHG,idigit, 0, 0, cellEFrac*cellAmplitude);
437 
438  if(nLabels > 0)
439  {
440  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
441  }
442 
443  if (fSubBackground)
444  {
445  Float_t energy = cellAmplitude;
446  Float_t time = cellTime;
447  fClusterizer->Calibrate(energy,time,cellNumber);
448  digit->SetAmplitude(energy);
449  avgE += energy;
450  }
451 
452  idigit++;
453  }
454 
455  if (fSubBackground) {
456  avgE /= fGeom->GetNumberOfSuperModules()*48*24;
457  Int_t ndigis = fDigitsArr->GetEntries();
458  for (Int_t i = 0; i < ndigis; ++i) {
459  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
460  Double_t energy = digit->GetAmplitude() - avgE;
461  if (energy<=0.001) {
462  digit->SetAmplitude(0);
463  } else {
464  digit->SetAmplitude(energy);
465  }
466  }
467  }
468  }
469 }
470 
476 {
477  const Int_t Ncls = fClusterArr->GetEntries();
478  AliDebug(1, Form("total no of clusters %d", Ncls));
479 
480  for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i)
481  {
482  AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
483 
484  Int_t ncellsTrue = 0;
485  const Int_t ncells = recpoint->GetMultiplicity();
486  UShort_t absIds[ncells];
487  Double32_t ratios[ncells];
488  Int_t *dlist = recpoint->GetDigitsList();
489  Float_t *elist = recpoint->GetEnergiesList();
490  Double_t mcEnergy = 0;
491 
492  for (Int_t c = 0; c < ncells; ++c)
493  {
494  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
495  absIds[ncellsTrue] = digit->GetId();
496  ratios[ncellsTrue] = elist[c]/digit->GetAmplitude();
497 
498  if (digit->GetIparent(1) > 0)
499  mcEnergy += digit->GetDEParent(1)/recpoint->GetEnergy();
500 
501  ++ncellsTrue;
502  }
503 
504  if (ncellsTrue < 1)
505  {
506  AliWarning("Skipping cluster with no cells");
507  continue;
508  }
509 
510  // calculate new cluster position
511  TVector3 gpos;
512  recpoint->GetGlobalPosition(gpos);
513  Float_t g[3];
514  gpos.GetXYZ(g);
515 
516  AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
517 
518  AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
519  c->SetType(AliVCluster::kEMCALClusterv1);
520  c->SetE(recpoint->GetEnergy());
521  c->SetPosition(g);
522  c->SetNCells(ncellsTrue);
523  c->SetCellsAbsId(absIds);
524  c->SetCellsAmplitudeFraction(ratios);
525  c->SetID(nout-1);
526  c->SetDispersion(recpoint->GetDispersion());
527  c->SetEmcCpvDistance(-1);
528  c->SetChi2(-1);
529  c->SetTOF(recpoint->GetTime()) ; //time-of-flight
530  c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
531  Float_t elipAxis[2];
532  recpoint->GetElipsAxis(elipAxis);
533  c->SetM02(elipAxis[0]*elipAxis[0]);
534  c->SetM20(elipAxis[1]*elipAxis[1]);
535  c->SetMCEnergyFraction(mcEnergy);
536 
537  //
538  // MC labels
539  //
540  Int_t parentMult = 0;
541  Int_t *parentList = recpoint->GetParents(parentMult);
542  Float_t *parentListDE = recpoint->GetParentsDE(); // deposited energy
543 
544  c->SetLabel(parentList, parentMult);
546  c->SetClusterMCEdepFractionFromEdepArray(parentListDE);
547  }
548 
549  //
550  // Set the cell energy deposition fraction map:
551  //
552  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
553  {
554  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
555 
556  // Get the digit that originated this cell cluster
557  //AliVCaloCells* cells = InputEvent()->GetEMCALCells();
558 
559  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
560  {
561  Int_t idigit = fCellLabels[absIds[icell]];
562 
563  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
564 
565  // Find the 4 MC labels that contributed to the cluster and their
566  // deposited energy in the current digit
567 
568  mcEdepFracPerCell[icell] = 0; // init
569 
570  Int_t nparents = dig->GetNiparent();
571  if ( nparents > 0 )
572  {
573  Int_t digLabel =-1 ;
574  Float_t edep = 0 ;
575  Float_t edepTot = 0 ;
576  Float_t mcEDepFrac[4] = {0,0,0,0};
577 
578  // all parents in digit
579  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
580  {
581  digLabel = dig->GetIparent (jndex+1);
582  edep = dig->GetDEParent(jndex+1);
583  edepTot += edep;
584 
585  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
586  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
587  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
588  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
589  } // all prarents in digit
590 
591  // Divide energy deposit by total deposited energy
592  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
593  if(edepTot > 0.01)
594  {
595  mcEdepFracPerCell[icell] = c->PackMCEdepFraction(mcEDepFrac);
596  }
597  } // at least one parent label in digit
598  } // cell in cluster loop
599 
600  c->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
601 
602  delete [] mcEdepFracPerCell;
603 
604  } // at least one parent in cluster, do the cell primary packing
605  }
606 }
607 
612 {
613  // Before destroying the orignal list, assign to the rec points the MC labels
614  // of the original clusters, if requested
617 
619 
620  fCaloClusters->Compress();
621 
623 }
624 
629 {
630  Int_t nclusters = fCaloClusters->GetEntriesFast();
631  for (Int_t icluster=0; icluster < nclusters; ++icluster) {
632  AliVCluster *clust = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
633  if (!clust) {
634  continue;
635  }
636  if (!clust->IsEMCAL()) {
637  continue;
638  }
639 
640  // SHOWER SHAPE -----------------------------------------------
641  if (fRecalShowerShape)
643 
644  // DISTANCE TO BAD CHANNELS -----------------------------------
647 
648  }
649 
650  fCaloClusters->Compress();
651 }
652 
657 {
658  if (label < 0) return;
659 
660  TClonesArray * arr = dynamic_cast<TClonesArray*>(fAod->FindListObject("mcparticles")) ;
661  if (!arr) return ;
662 
663  if (label < arr->GetEntriesFast())
664  {
665  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
666  if (!particle) return ;
667 
668  if (label == particle->Label()) return ; // label already OK
669  }
670 
671  // loop on the particles list and check if there is one with the same label
672  for (Int_t ind = 0; ind < arr->GetEntriesFast(); ind++)
673  {
674  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
675  if (!particle) continue ;
676 
677  if (label == particle->Label())
678  {
679  label = ind;
680  return;
681  }
682  }
683 
684  label = -1;
685 }
686 
696 {
697  Int_t ncls = fClusterArr->GetEntriesFast();
698  for(Int_t irp=0; irp < ncls; ++irp)
699  {
700  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
701  clArray.Reset();
702  Int_t nClu = 0;
703  Int_t nLabTotOrg = 0;
704  Float_t emax = -1;
705  Int_t idMax = -1;
706 
707  AliEMCALRecPoint *clus = static_cast<AliEMCALRecPoint*>(fClusterArr->At(irp));
708 
709  //Find the clusters that originally had the cells
710  const Int_t ncells = clus->GetMultiplicity();
711  Int_t *digList = clus->GetDigitsList();
712 
713  for (Int_t iLoopCell = 0 ; iLoopCell < ncells ; iLoopCell++)
714  {
715  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(digList[iLoopCell]));
716  Int_t idCell = digit->GetId();
717 
718  if (idCell>=0)
719  {
720  Int_t idCluster = fOrgClusterCellId[idCell];
721  Bool_t set = kTRUE;
722  for (Int_t icl =0; icl < nClu; icl++)
723  {
724  if (((Int_t)clArray.GetAt(icl))==-1) continue;
725  if (idCluster == ((Int_t)clArray.GetAt(icl))) set = kFALSE;
726  }
727  if (set && idCluster >= 0)
728  {
729  clArray.SetAt(idCluster,nClu++);
730  nLabTotOrg+=(fEventManager.InputEvent()->GetCaloCluster(idCluster))->GetNLabels();
731 
732  //Search highest E cluster
733  AliVCluster * clOrg = fEventManager.InputEvent()->GetCaloCluster(idCluster);
734  if (emax < clOrg->E())
735  {
736  emax = clOrg->E();
737  idMax = idCluster;
738  }
739  }
740  }
741  }// cell loop
742 
743  // Put the first in the list the cluster with highest energy
744  if (idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
745  {
746  Int_t maxIndex = -1;
747  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
748  for (Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++)
749  {
750  if (idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
751  }
752 
753  if (firstCluster >=0 && idMax >=0)
754  {
755  clArray.SetAt(idMax,0);
756  clArray.SetAt(firstCluster,maxIndex);
757  }
758  }
759 
760  // Get the labels list in the original clusters, assign all to the new cluster
761  TArrayI clMCArray(nLabTotOrg) ;
762  clMCArray.Reset();
763 
764  Int_t nLabTot = 0;
765  for (Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++)
766  {
767  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
768  AliVCluster * clOrg = fEventManager.InputEvent()->GetCaloCluster(idCluster);
769  Int_t nLab = clOrg->GetNLabels();
770 
771  for (Int_t iLab = 0 ; iLab < nLab ; iLab++)
772  {
773  Int_t lab = clOrg->GetLabelAt(iLab) ;
774  if (lab>=0)
775  {
776  Bool_t set = kTRUE;
777  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
778  {
779  if (lab == ((Int_t)clMCArray.GetAt(iLabTot))) set = kFALSE;
780  }
781  if (set) clMCArray.SetAt(lab,nLabTot++);
782  }
783  }
784  }// cluster loop
785 
786  // Set the final list of labels to rec point
787 
788  Int_t *labels = new Int_t[nLabTot];
789  for(Int_t il = 0; il < nLabTot; il++) labels[il] = clMCArray.GetArray()[il];
790  clus->SetParents(nLabTot,labels);
791 
792  } // rec point array
793 }
794 
799 {
800  // Select clusterization/unfolding algorithm and set all the needed parameters.
801 
802  // If distBC option enabled, init and fill bad channel map
805  else
807 
809 
810  if (fJustUnfold){
811  // init the unfolding afterburner
812  delete fUnfolder;
813  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
814  return;
815  }
816 
817  if (!fGeomMatrixSet) {
818  if (fLoadGeomMatrices) {
819  for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
820  if (fGeomMatrix[mod]){
821  if (AliAnalysisManager::GetAnalysisManager()->GetDebugLevel() > 2)
822  fGeomMatrix[mod]->Print();
823  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
824  }
825  }
826  } else { // get matrix from file (work around bug in aliroot)
827  for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
828  const TGeoHMatrix *gm = 0;
829  if (fEsd) {
830  gm = fEsd->GetEMCALMatrix(mod);
831  } else {
832  AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(fAod->GetHeader());
833  if(!aodheader) AliFatal("Not a standard AOD");
834  if (aodheader) {
835  gm = aodheader->GetEMCALMatrix(mod);
836  }
837  }
838  if (gm) {
839  if (AliAnalysisManager::GetAnalysisManager()->GetDebugLevel() > 2)
840  gm->Print();
841  fGeom->SetMisalMatrix(gm,mod);
842  }
843  }
844  }
845  fGeomMatrixSet=kTRUE;
846  }
847 
848  // setup digit array if needed
849  if (!fDigitsArr) {
850  fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
851  fDigitsArr->SetOwner(1);
852  }
853 
854  // then setup clusterizer
855  if (fClusterizer) {
856  // avoid to delete digits array
857  fClusterizer->SetDigitsArr(0);
858  delete fClusterizer;
859  }
860  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
861  fClusterizer = new AliEMCALClusterizerv1(fGeom);
862  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
863  AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
864  clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
865  clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
866  fClusterizer = clusterizer;
867  }
868  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
869  fClusterizer = new AliEMCALClusterizerv2(fGeom);
870  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv3)
871  fClusterizer = new AliEMCALClusterizerv3(fGeom);
872  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
873  AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
874  clusterizer->SetNphi(fNPhi);
875  clusterizer->SetNeta(fNEta);
876  clusterizer->SetShiftPhi(fShiftPhi);
877  clusterizer->SetShiftEta(fShiftEta);
878  clusterizer->SetTRUshift(fTRUShift);
879  fClusterizer = clusterizer;
880  }
881  else {
882  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
883  }
884  fClusterizer->InitParameters(fRecParam);
885 
886  if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
887  AliCDBManager *cdb = AliCDBManager::Instance();
888  if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
889  cdb->SetDefaultStorage(fOCDBpath);
890  if (fRun!=cdb->GetRun())
891  cdb->SetRun(fRun);
892  }
893  if (!fCalibData&&fLoadCalib&&fRun>0) {
894  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
895  if (entry)
896  fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
897  if (!fCalibData)
898  AliFatal("Calibration parameters not found in CDB!");
899  }
900  if (!fPedestalData&&fLoadPed&&fRun>0) {
901  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
902  if (entry)
903  fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
904  }
905  if (fCalibData) {
906  fClusterizer->SetInputCalibrated(kFALSE);
907  fClusterizer->SetCalibrationParameters(fCalibData);
908  } else {
909  fClusterizer->SetInputCalibrated(kTRUE);
910  }
911  fClusterizer->SetCaloCalibPedestal(fPedestalData);
912  fClusterizer->SetJustClusters(kTRUE);
913  fClusterizer->SetDigitsArr(fDigitsArr);
914  fClusterizer->SetOutput(0);
915  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
916 
917 }
918 
924 {
926 
927  if (runChanged && fRecalDistToBadChannels) {
928  // init bad channels
929  Int_t fInitBC = InitBadChannels();
930  if (fInitBC==0) {
931  AliError("InitBadChannels returned false, returning");
932  }
933  if (fInitBC==1) {
934  AliWarning("InitBadChannels OK");
935  }
936  if (fInitBC>1) {
937  AliWarning(Form("No external hot channel set: %d - %s", fEventManager.InputEvent()->GetRunNumber(), fFilepass.Data()));
938  }
939  }
940  return runChanged;
941 }
942 
947 {
948  const Int_t nents = fCaloClusters->GetEntries();
949  for (Int_t i=0;i<nents;++i) {
950  AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
951  if (!c) {
952  continue;
953  }
954  if (c->IsEMCAL()) {
955  delete fCaloClusters->RemoveAt(i);
956  }
957  }
958 }
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:44
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
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
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.
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.
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
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)