AliPhysics  63e47e1 (63e47e1)
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  return kTRUE;
210 }
211 
216 {
218 
219  if (fCreateHisto){
220  fHistCPUTime = new TH1F("hCPUTime","hCPUTime;CPU Time (ms)", 2000, 0, 1000);
221  fOutput->Add(fHistCPUTime);
222 
223  fHistRealTime = new TH1F("hRealTime","hRealTime;Real Time (ms)", 2000, 0, 1000);
224  fOutput->Add(fHistRealTime);
225 
226  fTimer = new TStopwatch();
227  }
228 }
229 
234 {
235  // Time the event loop if histogram creation is enabled
236  if (fCreateHisto) {
237  fTimer->Start(kTRUE);
238  }
239 
241 
242  fEsd = dynamic_cast<AliESDEvent*>(fEventManager.InputEvent());
243  fAod = dynamic_cast<AliAODEvent*>(fEventManager.InputEvent());
244 
245  // Only support one cluster container in the clusterizer!
247  if (!clusCont) {
248  AliFatal("Could not retrieve cluster container!");
249  }
250 
251  fCaloClusters = clusCont->GetArray();
252 
253  // If cells are empty, clear clusters and return
254  if (fCaloCells->GetNumberOfCells()<=0)
255  {
256  AliWarning(Form("Number of EMCAL cells = %d, returning", fCaloCells->GetNumberOfCells()));
258  return kFALSE;
259  }
260 
261  UInt_t offtrigger = 0;
262  if (fEsd) {
263  UInt_t mask1 = fEsd->GetESDRun()->GetDetectorsInDAQ();
264  UInt_t mask2 = fEsd->GetESDRun()->GetDetectorsInReco();
265  Bool_t desc1 = (mask1 >> 18) & 0x1;
266  Bool_t desc2 = (mask2 >> 18) & 0x1;
267  if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
268  AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
269  mask1, fEsd->GetESDRun()->GetDetectorsInReco(),
270  mask2, fEsd->GetESDRun()->GetDetectorsInDAQ()));
271  return kFALSE;
272  }
273  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
274  offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
275  } else {
276  offtrigger = ((AliVAODHeader*)fAod->GetHeader())->GetOfflineTrigger();
277  }
278 
279  if (!fMCEvent) {
280  if (offtrigger & AliVEvent::kFastOnly) {
281  AliError(Form("EMCAL not in fast only partition"));
282  return kFALSE;
283  }
284  }
285 
286  Init();
287 
288  if (fJustUnfold) {
289  AliWarning("Unfolding not implemented");
290  return kTRUE;
291  }
292 
293  FillDigitsArray();
294 
295  Clusterize();
296 
297  UpdateClusters();
298 
300 
301  if (fCreateHisto) {
302  fTimer->Stop();
303  // Ensure that it is stored in ms
304  fHistCPUTime->Fill(fTimer->CpuTime() * 1000.);
305  fHistRealTime->Fill(fTimer->RealTime() * 1000.);
306  }
307 
308  return kTRUE;
309 }
310 
315 {
316  if (fSubBackground) {
317  fClusterizer->SetInputCalibrated(kTRUE);
318  fClusterizer->SetCalibrationParameters(0);
319  }
320 
321  fClusterizer->Digits2Clusters("");
322 
323  if (fSubBackground) {
324  if (fCalibData) {
325  fClusterizer->SetInputCalibrated(kFALSE);
326  fClusterizer->SetCalibrationParameters(fCalibData);
327  }
328  }
329 }
330 
335 {
336  fDigitsArr->Clear("C");
337  if (fTestPatternInput) {
338  // Fill digits from a pattern
339  Int_t maxd = fGeom->GetNCells() / 4;
340  for (Int_t idigit = 0; idigit < maxd; idigit++){
341  if (idigit % 24 == 12) idigit += 12;
342  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
343  digit->SetId(idigit * 4);
344  digit->SetTime(600);
345  digit->SetTimeR(600);
346  digit->SetIndexInList(idigit);
347  digit->SetType(AliEMCALDigit::kHG);
348  digit->SetAmplitude(0.1);
349  }
350  }
351  else
352  {
353  // In case of MC productions done before aliroot tag v5-02-Rev09
354  // passing the cluster label to all the cells belonging to this cluster
355  // very rough
356  // Copied and simplified from AliEMCALTenderSupply
358  {
359  for (Int_t i = 0; i < fgkTotalCellNumber; i++)
360  {
361  fCellLabels [i] =-1 ;
362  fOrgClusterCellId[i] =-1 ;
363  }
364 
365  Int_t nClusters = fEventManager.InputEvent()->GetNumberOfCaloClusters();
366  for (Int_t i = 0; i < nClusters; i++)
367  {
368  AliVCluster *clus = fEventManager.InputEvent()->GetCaloCluster(i);
369 
370  if (!clus) continue;
371 
372  if (!clus->IsEMCAL()) continue ;
373 
374  Int_t label = clus->GetLabel();
375  UShort_t * index = clus->GetCellsAbsId() ;
376 
377  for(Int_t icell=0; icell < clus->GetNCells(); icell++)
378  {
380  fCellLabels[index[icell]] = label;
381 
382  fOrgClusterCellId[index[icell]] = i ; // index of the original cluster
383  } // cell in cluster loop
384  } // cluster loop
385  }
386 
387  Double_t avgE = 0; // for background subtraction
388  const Int_t ncells = fCaloCells->GetNumberOfCells();
389  for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell)
390  {
391  Float_t cellAmplitude=0;
392  Double_t cellTime=0, amp = 0, cellEFrac = 0;
393  Short_t cellNumber=0;
394  Int_t cellMCLabel=-1;
395  if (fCaloCells->GetCell(icell, cellNumber, amp, cellTime, cellMCLabel, cellEFrac) != kTRUE)
396  break;
397 
398  cellAmplitude = amp; // compilation problem
399 
400  //if (fSetCellMCLabelFromCluster) cellMCLabel = fCellLabels[cellNumber];
402  {
403  if (fSetCellMCLabelFromCluster) cellMCLabel = fCellLabels[cellNumber];
404  else if (fRemapMCLabelForAODs ) RemapMCLabelForAODs(cellMCLabel);
405  }
406 
407  if (cellMCLabel > 0 && cellEFrac < 1e-6)
408  cellEFrac = 1;
409 
410  if (cellAmplitude < 1e-6 || cellNumber < 0)
411  continue;
412 
413  // New way to set the cell MC labels,
414  // valid only for MC productions with aliroot > v5-07-21
415  //
416  TArrayI labeArr(0);
417  TArrayF eDepArr(0);
418  Int_t nLabels = 0;
419 
420  if(fSetCellMCLabelFromEdepFrac && fOrgClusterCellId[cellNumber] >= 0) // index can be negative if noisy cell that did not form cluster
421  {
422  cellMCLabel = -1;
423 
424  fCellLabels[cellNumber] = idigit;
425 
426  Int_t iclus = fOrgClusterCellId[cellNumber];
427 
428  if(iclus < 0)
429  {
430  AliInfo("Negative original cluster index, skip \n");
431  continue;
432  }
433 
434  AliVCluster* clus = fEventManager.InputEvent()->GetCaloCluster(iclus);
435 
436  for(Int_t icluscell=0; icluscell < clus->GetNCells(); icluscell++ )
437  {
438  if(cellNumber != clus->GetCellAbsId(icluscell)) continue ;
439 
440  // Select the MC label contributing, only if enough energy deposition
441  fRecoUtils->RecalculateCellLabelsRemoveAddedGenerator(cellNumber, clus, fMCEvent, cellAmplitude, labeArr, eDepArr);
442  nLabels = labeArr.GetSize();
443  }
444  }
445 
446  AliEMCALDigit *digit = new((*fDigitsArr)[idigit]) AliEMCALDigit(cellMCLabel, cellMCLabel, cellNumber,
447  cellAmplitude, (Float_t)cellTime,
448  AliEMCALDigit::kHG,idigit, 0, 0, cellEFrac*cellAmplitude);
449 
450  if(nLabels > 0)
451  {
452  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
453  }
454 
455  if (fSubBackground)
456  {
457  Float_t energy = cellAmplitude;
458  Float_t time = cellTime;
459  fClusterizer->Calibrate(energy,time,cellNumber);
460  digit->SetAmplitude(energy);
461  avgE += energy;
462  }
463 
464  idigit++;
465  }
466 
467  if (fSubBackground) {
468  avgE /= fGeom->GetNumberOfSuperModules()*48*24;
469  Int_t ndigis = fDigitsArr->GetEntries();
470  for (Int_t i = 0; i < ndigis; ++i) {
471  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
472  Double_t energy = digit->GetAmplitude() - avgE;
473  if (energy<=0.001) {
474  digit->SetAmplitude(0);
475  } else {
476  digit->SetAmplitude(energy);
477  }
478  }
479  }
480  }
481 }
482 
488 {
489  const Int_t Ncls = fClusterArr->GetEntries();
490  AliDebug(1, Form("total no of clusters %d", Ncls));
491 
492  for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i)
493  {
494  AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
495 
496  Int_t ncellsTrue = 0;
497  const Int_t ncells = recpoint->GetMultiplicity();
498  UShort_t absIds[ncells];
499  Double32_t ratios[ncells];
500  Int_t *dlist = recpoint->GetDigitsList();
501  Float_t *elist = recpoint->GetEnergiesList();
502  Double_t mcEnergy = 0;
503 
504  Float_t clusterE = 0;
505  for (Int_t c = 0; c < ncells; ++c)
506  {
507  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
508  absIds[ncellsTrue] = digit->GetId();
509  ratios[ncellsTrue] = elist[c]/digit->GetAmplitude();
510 
511  if ( !fRecParam->GetUnfold() && (ratios[ncellsTrue] > 1 || ratios[ncellsTrue] < 1) )
512  AliWarning(Form("recpoint cell E %2.3f but digit E %2.3f and no unfolding",
513  recpoint->GetEnergiesList()[c], digit->GetAmplitude()));
514 
515  // In case of unfolding, remove digits with unfolded energy too low
516  if ( fRecParam->GetUnfold() )
517  {
518  if ( recpoint->GetEnergiesList()[c] < fUnfoldCellMinE ||
519  ratios[ncellsTrue] < fUnfoldCellMinEFrac )
520  {
521 
522  AliDebug(2,Form("Too small energy in cell of cluster: cluster cell %f, digit %f",
523  recpoint->GetEnergiesList()[c],digit->GetAmplitude()));
524  continue;
525  } // if cuts
526  }// Select cells
527 
528  // Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
529  clusterE += recpoint->GetEnergiesList()[c];
530 
531  if (digit->GetIparent(1) > 0)
532  mcEnergy += digit->GetDEParent(1)/recpoint->GetEnergy();
533 
534  ++ncellsTrue;
535  } // cluster cell loop
536 
537  if (ncellsTrue < 1)
538  {
539  AliWarning("Skipping cluster with no cells");
540  continue;
541  }
542 
543  // calculate new cluster position
544  TVector3 gpos;
545  recpoint->GetGlobalPosition(gpos);
546  Float_t g[3];
547  gpos.GetXYZ(g);
548 
549  AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
550 
551  AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
552  c->SetType(AliVCluster::kEMCALClusterv1);
553  c->SetE(clusterE);
554  c->SetPosition(g);
555  c->SetNCells(ncellsTrue);
556  c->SetCellsAbsId(absIds);
557  c->SetCellsAmplitudeFraction(ratios);
558  c->SetID(nout-1);
559  c->SetDispersion(recpoint->GetDispersion());
560  c->SetEmcCpvDistance(-1);
561  c->SetChi2(-1);
562  c->SetTOF(recpoint->GetTime()) ; //time-of-flight
563  c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
564  c->SetMCEnergyFraction(mcEnergy);
565 
566  if ( ncells == ncellsTrue )
567  {
568  Float_t elipAxis[2];
569  recpoint->GetElipsAxis(elipAxis);
570  c->SetM02(elipAxis[0]*elipAxis[0]);
571  c->SetM20(elipAxis[1]*elipAxis[1]);
572  }
573  else
574  {
575  // In case some cells rejected, in unfolding case, recalculate
576  // shower shape parameters and position
577  AliDebug(2,Form("Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape",ncells,ncellsTrue));
578 
580  //fRecoUtils->RecalculateClusterPID(c);
582  }
583 
584  //
585  // MC labels
586  //
587  Int_t parentMult = 0;
588  Int_t *parentList = recpoint->GetParents(parentMult);
589  Float_t *parentListDE = recpoint->GetParentsDE(); // deposited energy
590 
591  c->SetLabel(parentList, parentMult);
593  c->SetClusterMCEdepFractionFromEdepArray(parentListDE);
594  }
595 
596  //
597  // Set the cell energy deposition fraction map:
598  //
599  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
600  {
601  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
602 
603  // Get the digit that originated this cell cluster
604  //AliVCaloCells* cells = InputEvent()->GetEMCALCells();
605 
606  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
607  {
608  Int_t idigit = fCellLabels[absIds[icell]];
609 
610  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
611 
612  // Find the 4 MC labels that contributed to the cluster and their
613  // deposited energy in the current digit
614 
615  mcEdepFracPerCell[icell] = 0; // init
616 
617  Int_t nparents = dig->GetNiparent();
618  if ( nparents > 0 )
619  {
620  Int_t digLabel =-1 ;
621  Float_t edep = 0 ;
622  Float_t edepTot = 0 ;
623  Float_t mcEDepFrac[4] = {0,0,0,0};
624 
625  // all parents in digit
626  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
627  {
628  digLabel = dig->GetIparent (jndex+1);
629  edep = dig->GetDEParent(jndex+1);
630  edepTot += edep;
631 
632  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
633  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
634  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
635  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
636  } // all prarents in digit
637 
638  // Divide energy deposit by total deposited energy
639  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
640  if(edepTot > 0.01)
641  {
642  mcEdepFracPerCell[icell] = c->PackMCEdepFraction(mcEDepFrac);
643  }
644  } // at least one parent label in digit
645  } // cell in cluster loop
646 
647  c->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
648 
649  delete [] mcEdepFracPerCell;
650 
651  } // at least one parent in cluster, do the cell primary packing
652  }
653 }
654 
659 {
660  // Before destroying the orignal list, assign to the rec points the MC labels
661  // of the original clusters, if requested
664 
666 
667  fCaloClusters->Compress();
668 
670 }
671 
676 {
677  Int_t nclusters = fCaloClusters->GetEntriesFast();
678  for (Int_t icluster=0; icluster < nclusters; ++icluster) {
679  AliVCluster *clust = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
680  if (!clust) {
681  continue;
682  }
683  if (!clust->IsEMCAL()) {
684  continue;
685  }
686 
687  // SHOWER SHAPE -----------------------------------------------
688  if (fRecalShowerShape)
690 
691  // DISTANCE TO BAD CHANNELS -----------------------------------
694 
695  }
696 
697  fCaloClusters->Compress();
698 }
699 
704 {
705  if (label < 0) return;
706 
707  TClonesArray * arr = dynamic_cast<TClonesArray*>(fAod->FindListObject("mcparticles")) ;
708  if (!arr) return ;
709 
710  if (label < arr->GetEntriesFast())
711  {
712  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
713  if (!particle) return ;
714 
715  if (label == particle->Label()) return ; // label already OK
716  }
717 
718  // loop on the particles list and check if there is one with the same label
719  for (Int_t ind = 0; ind < arr->GetEntriesFast(); ind++)
720  {
721  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
722  if (!particle) continue ;
723 
724  if (label == particle->Label())
725  {
726  label = ind;
727  return;
728  }
729  }
730 
731  label = -1;
732 }
733 
743 {
744  Int_t ncls = fClusterArr->GetEntriesFast();
745  for(Int_t irp=0; irp < ncls; ++irp)
746  {
747  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
748  clArray.Reset();
749  Int_t nClu = 0;
750  Int_t nLabTotOrg = 0;
751  Float_t emax = -1;
752  Int_t idMax = -1;
753 
754  AliEMCALRecPoint *clus = static_cast<AliEMCALRecPoint*>(fClusterArr->At(irp));
755 
756  //Find the clusters that originally had the cells
757  const Int_t ncells = clus->GetMultiplicity();
758  Int_t *digList = clus->GetDigitsList();
759 
760  for (Int_t iLoopCell = 0 ; iLoopCell < ncells ; iLoopCell++)
761  {
762  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(digList[iLoopCell]));
763  Int_t idCell = digit->GetId();
764 
765  if (idCell>=0)
766  {
767  Int_t idCluster = fOrgClusterCellId[idCell];
768  Bool_t set = kTRUE;
769  for (Int_t icl =0; icl < nClu; icl++)
770  {
771  if (((Int_t)clArray.GetAt(icl))==-1) continue;
772  if (idCluster == ((Int_t)clArray.GetAt(icl))) set = kFALSE;
773  }
774  if (set && idCluster >= 0)
775  {
776  clArray.SetAt(idCluster,nClu++);
777  nLabTotOrg+=(fEventManager.InputEvent()->GetCaloCluster(idCluster))->GetNLabels();
778 
779  //Search highest E cluster
780  AliVCluster * clOrg = fEventManager.InputEvent()->GetCaloCluster(idCluster);
781  if (emax < clOrg->E())
782  {
783  emax = clOrg->E();
784  idMax = idCluster;
785  }
786  }
787  }
788  }// cell loop
789 
790  // Put the first in the list the cluster with highest energy
791  if (idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
792  {
793  Int_t maxIndex = -1;
794  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
795  for (Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++)
796  {
797  if (idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
798  }
799 
800  if (firstCluster >=0 && idMax >=0)
801  {
802  clArray.SetAt(idMax,0);
803  clArray.SetAt(firstCluster,maxIndex);
804  }
805  }
806 
807  // Get the labels list in the original clusters, assign all to the new cluster
808  TArrayI clMCArray(nLabTotOrg) ;
809  clMCArray.Reset();
810 
811  Int_t nLabTot = 0;
812  for (Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++)
813  {
814  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
815  AliVCluster * clOrg = fEventManager.InputEvent()->GetCaloCluster(idCluster);
816  Int_t nLab = clOrg->GetNLabels();
817 
818  for (Int_t iLab = 0 ; iLab < nLab ; iLab++)
819  {
820  Int_t lab = clOrg->GetLabelAt(iLab) ;
821  if (lab>=0)
822  {
823  Bool_t set = kTRUE;
824  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
825  {
826  if (lab == ((Int_t)clMCArray.GetAt(iLabTot))) set = kFALSE;
827  }
828  if (set) clMCArray.SetAt(lab,nLabTot++);
829  }
830  }
831  }// cluster loop
832 
833  // Set the final list of labels to rec point
834 
835  Int_t *labels = new Int_t[nLabTot];
836  for(Int_t il = 0; il < nLabTot; il++) labels[il] = clMCArray.GetArray()[il];
837  clus->SetParents(nLabTot,labels);
838 
839  } // rec point array
840 }
841 
846 {
847  // Select clusterization/unfolding algorithm and set all the needed parameters.
848 
849  // If distBC option enabled, init and fill bad channel map
852  else
854 
856 
857  if (fJustUnfold){
858  // init the unfolding afterburner
859  delete fUnfolder;
860  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
861  return;
862  }
863 
864  if (!fGeomMatrixSet) {
865  if (fLoadGeomMatrices) {
866  for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
867  if (fGeomMatrix[mod]){
868  if (AliAnalysisManager::GetAnalysisManager()->GetDebugLevel() > 2)
869  fGeomMatrix[mod]->Print();
870  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
871  }
872  }
873  } else { // get matrix from file (work around bug in aliroot)
874  for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
875  const TGeoHMatrix *gm = 0;
876  if (fEsd) {
877  gm = fEsd->GetEMCALMatrix(mod);
878  } else {
879  AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(fAod->GetHeader());
880  if(!aodheader) AliFatal("Not a standard AOD");
881  if (aodheader) {
882  gm = aodheader->GetEMCALMatrix(mod);
883  }
884  }
885  if (gm) {
886  if (AliAnalysisManager::GetAnalysisManager()->GetDebugLevel() > 2)
887  gm->Print();
888  fGeom->SetMisalMatrix(gm,mod);
889  }
890  }
891  }
892  fGeomMatrixSet=kTRUE;
893  }
894 
895  // setup digit array if needed
896  if (!fDigitsArr) {
897  fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
898  fDigitsArr->SetOwner(1);
899  }
900 
901  // then setup clusterizer
902  if (fClusterizer) {
903  // avoid to delete digits array
904  fClusterizer->SetDigitsArr(0);
905  delete fClusterizer;
906  }
907  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
908  fClusterizer = new AliEMCALClusterizerv1(fGeom);
909  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
910  AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
911  clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
912  clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
913  fClusterizer = clusterizer;
914  }
915  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
916  fClusterizer = new AliEMCALClusterizerv2(fGeom);
917  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv3)
918  fClusterizer = new AliEMCALClusterizerv3(fGeom);
919  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
920  AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
921  clusterizer->SetNphi(fNPhi);
922  clusterizer->SetNeta(fNEta);
923  clusterizer->SetShiftPhi(fShiftPhi);
924  clusterizer->SetShiftEta(fShiftEta);
925  clusterizer->SetTRUshift(fTRUShift);
926  fClusterizer = clusterizer;
927  }
928  else {
929  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
930  }
931  fClusterizer->InitParameters(fRecParam);
932 
933  if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
934  AliCDBManager *cdb = AliCDBManager::Instance();
935  if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
936  cdb->SetDefaultStorage(fOCDBpath);
937  if (fRun!=cdb->GetRun())
938  cdb->SetRun(fRun);
939  }
940  if (!fCalibData&&fLoadCalib&&fRun>0) {
941  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
942  if (entry)
943  fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
944  if (!fCalibData)
945  AliFatal("Calibration parameters not found in CDB!");
946  }
947  if (!fPedestalData&&fLoadPed&&fRun>0) {
948  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
949  if (entry)
950  fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
951  }
952  if (fCalibData) {
953  fClusterizer->SetInputCalibrated(kFALSE);
954  fClusterizer->SetCalibrationParameters(fCalibData);
955  } else {
956  fClusterizer->SetInputCalibrated(kTRUE);
957  }
958  fClusterizer->SetCaloCalibPedestal(fPedestalData);
959  fClusterizer->SetJustClusters(kTRUE);
960  fClusterizer->SetDigitsArr(fDigitsArr);
961  fClusterizer->SetOutput(0);
962  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
963 
964  //fClusterizer->Print("");
965 }
966 
972 {
974 
975  if (runChanged && fRecalDistToBadChannels) {
976  // init bad channels
977  Int_t fInitBC = InitBadChannels();
978  if (fInitBC==0) {
979  AliError("InitBadChannels returned false, returning");
980  }
981  if (fInitBC==1) {
982  AliWarning("InitBadChannels OK");
983  }
984  if (fInitBC>1) {
985  AliWarning(Form("No external hot channel set: %d - %s", fEventManager.InputEvent()->GetRunNumber(), fFilepass.Data()));
986  }
987  }
988  return runChanged;
989 }
990 
995 {
996  const Int_t nents = fCaloClusters->GetEntries();
997  for (Int_t i=0;i<nents;++i) {
998  AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
999  if (!c) {
1000  continue;
1001  }
1002  if (c->IsEMCAL()) {
1003  delete fCaloClusters->RemoveAt(i);
1004  }
1005  }
1006 }
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
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.
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
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)