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