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