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