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