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