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