AliPhysics  59e0e03 (59e0e03)
AliAnalysisTaskEMCALClusterizeFast.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- Root ---
17 #include <TClonesArray.h>
18 #include <TGeoManager.h>
19 #include <TObjArray.h>
20 #include <TString.h>
21 #include <TTree.h>
22 #include <TArrayI.h>
23 
24 // --- AliRoot ---
25 #include "AliAODCaloCluster.h"
26 #include "AliAODEvent.h"
27 #include "AliAnalysisManager.h"
28 #include "AliCDBEntry.h"
29 #include "AliCDBManager.h"
30 #include "AliCaloCalibPedestal.h"
31 #include "AliEMCALAfterBurnerUF.h"
32 #include "AliEMCALCalibData.h"
33 #include "AliEMCALClusterizerNxN.h"
34 #include "AliEMCALClusterizerv1.h"
35 #include "AliEMCALClusterizerv2.h"
36 #include "AliEMCALClusterizerv3.h"
37 #include "AliEMCALClusterizerFixedWindow.h"
38 #include "AliEMCALDigit.h"
39 #include "AliEMCALGeometry.h"
40 #include "AliEMCALRecParam.h"
41 #include "AliEMCALRecPoint.h"
42 #include "AliEMCALRecoUtils.h"
43 #include "AliESDEvent.h"
44 #include "AliInputEventHandler.h"
45 #include "AliLog.h"
46 
48 
50 
51 //________________________________________________________________________
54  fRun(-1),
55  fDigitsArr(0),
56  fClusterArr(0),
57  fRecParam(new AliEMCALRecParam),
58  fClusterizer(0),
59  fUnfolder(0),
60  fJustUnfold(kFALSE),
61  fGeomName(),
62  fGeomMatrixSet(kFALSE),
63  fLoadGeomMatrices(kFALSE),
64  fOCDBpath(),
65  fCalibData(0),
66  fPedestalData(0),
67  fOutputAODBranch(0),
68  fOutputAODBrName(),
69  fRecoUtils(0),
70  fLoadCalib(kFALSE),
71  fLoadPed(kFALSE),
72  fAttachClusters(kTRUE),
73  fSubBackground(kFALSE),
74  fNPhi(4),
75  fNEta(4),
76  fShiftPhi(2),
77  fShiftEta(2),
78  fTRUShift(0),
79  fInputCellType(kFEEData),
80  fTrackName(),
81  fCaloCellsName(),
82  fCaloClustersName("newCaloClusters"),
83  fDoUpdateCells(kFALSE),
84  fDoClusterize(kTRUE),
85  fClusterBadChannelCheck(kFALSE),
86  fRejectExoticClusters(kFALSE),
87  fRejectExoticCells(kFALSE),
88  fFiducial(kFALSE),
89  fDoNonLinearity(kFALSE),
90  fRecalDistToBadChannels(kFALSE),
91  fSetCellMCLabelFromCluster(0),
92  fSetCellMCLabelFromEdepFrac(0),
93  fCaloCells(0),
94  fCaloClusters(0),
95  fEsd(0),
96  fAod(0),
97  fGeom(0)
98 {
99  // Constructor
100 
101  for(Int_t i = 0; i < AliEMCALGeoParams::fgkEMCALModules; i++) fGeomMatrix[i] = 0 ;
102  for(Int_t j = 0; j < fgkTotalCellNumber; j++)
103  { fOrgClusterCellId[j] =-1; fCellLabels[j] =-1 ; }
104 }
105 
106 //________________________________________________________________________
108  AliAnalysisTaskSE(name),
109  fRun(-1),
110  fDigitsArr(0),
111  fClusterArr(0),
112  fRecParam(new AliEMCALRecParam),
113  fClusterizer(0),
114  fUnfolder(0),
115  fJustUnfold(kFALSE),
116  fGeomName(),
117  fGeomMatrixSet(kFALSE),
118  fLoadGeomMatrices(kFALSE),
119  fOCDBpath(),
120  fCalibData(0),
121  fPedestalData(0),
122  fOutputAODBranch(0),
123  fOutputAODBrName(),
124  fRecoUtils(0),
125  fLoadCalib(kFALSE),
126  fLoadPed(kFALSE),
127  fAttachClusters(kTRUE),
128  fSubBackground(kFALSE),
129  fNPhi(4),
130  fNEta(4),
131  fShiftPhi(2),
132  fShiftEta(2),
133  fTRUShift(0),
134  fInputCellType(kFEEData),
135  fTrackName(),
136  fCaloCellsName(),
137  fCaloClustersName("newCaloClusters"),
138  fDoUpdateCells(kFALSE),
139  fDoClusterize(kTRUE),
140  fClusterBadChannelCheck(kFALSE),
141  fRejectExoticClusters(kFALSE),
142  fRejectExoticCells(kFALSE),
143  fFiducial(kFALSE),
144  fDoNonLinearity(kFALSE),
145  fRecalDistToBadChannels(kFALSE),
146  fSetCellMCLabelFromCluster(0),
147  fSetCellMCLabelFromEdepFrac(0),
148  fCaloCells(0),
149  fCaloClusters(0),
150  fEsd(0),
151  fAod(0),
152  fGeom(0)
153 {
154  // Constructor
155 
156  fBranchNames = "ESD:AliESDHeader.,AliESDRun.,EMCALCells.,EMCALTrigger. AOD:header,emcalCells";
157 
158  for(Int_t i = 0; i < AliEMCALGeoParams::fgkEMCALModules; i++) fGeomMatrix[i] = 0 ;
159  for(Int_t j = 0; j < fgkTotalCellNumber; j++)
160  { fOrgClusterCellId[j] =-1; fCellLabels[j] =-1 ; }
161 }
162 
163 //________________________________________________________________________
165 {
166  // Destructor.
167 
168  delete fClusterizer;
169  delete fUnfolder;
170  delete fRecoUtils;
171  delete fRecParam;
172 }
173 
174 //________________________________________________________________________
176 {
177  // Create output objects.
178 
179  if (!fOutputAODBrName.IsNull()) {
180  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
182  AddAODBranch("TClonesArray", &fOutputAODBranch);
183  AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
184  }
185 }
186 
187 //________________________________________________________________________
189 {
190  // Main loop, called for each event
191 
192  // if no operation is requested, return
193  if (!fDoClusterize && !fDoUpdateCells)
194  return;
195 
196  // remove the contents of output list set in the previous event
197  if (fOutputAODBranch)
198  fOutputAODBranch->Clear("C");
199 
200  fEsd = dynamic_cast<AliESDEvent*>(InputEvent());
201  fAod = dynamic_cast<AliAODEvent*>(InputEvent());
202 
203  if (!fEsd&&!fAod) {
204  Error("UserExec","Event not available");
205  return;
206  }
207 
208  LoadBranches();
209 
210  UInt_t offtrigger = 0;
211  if (fEsd) {
212  UInt_t mask1 = fEsd->GetESDRun()->GetDetectorsInDAQ();
213  UInt_t mask2 = fEsd->GetESDRun()->GetDetectorsInReco();
214  Bool_t desc1 = (mask1 >> 18) & 0x1;
215  Bool_t desc2 = (mask2 >> 18) & 0x1;
216  if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
217  AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
218  mask1, fEsd->GetESDRun()->GetDetectorsInReco(),
219  mask2, fEsd->GetESDRun()->GetDetectorsInDAQ()));
220  return;
221  }
222  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
223  offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
224  } else {
225  offtrigger = ((AliVAODHeader*)fAod->GetHeader())->GetOfflineTrigger();
226  }
227 
228  if (!MCEvent()) {
229  if (offtrigger & AliVEvent::kFastOnly) {
230  AliError(Form("EMCAL not in fast only partition"));
231  return;
232  }
233  }
234 
235  Init();
236 
237  if (fJustUnfold) {
238  AliWarning("Unfolding not implemented");
239  return;
240  }
241 
242  FillDigitsArray();
243 
244  if (fDoClusterize)
245  Clusterize();
246 
247  if (fDoUpdateCells)
248  UpdateCells();
249 
251  return;
252 
253  UpdateClusters();
255 
258 
259 }
260 
261 //________________________________________________________________________
262 void AliAnalysisTaskEMCALClusterizeFast::CopyClusters(TClonesArray *orig, TClonesArray *dest)
263 {
264  const Int_t Ncls = orig->GetEntries();
265 
266  for(Int_t i=0; i < Ncls; ++i) {
267  AliVCluster *oc = static_cast<AliVCluster*>(orig->At(i));
268 
269  if (!oc)
270  continue;
271 
272  if (!oc->IsEMCAL())
273  continue;
274 
275  AliVCluster *dc = static_cast<AliVCluster*>(dest->New(i));
276  dc->SetType(AliVCluster::kEMCALClusterv1);
277  dc->SetE(oc->E());
278  Float_t pos[3] = {0};
279  oc->GetPosition(pos);
280  dc->SetPosition(pos);
281  dc->SetNCells(oc->GetNCells());
282  dc->SetCellsAbsId(oc->GetCellsAbsId());
283  dc->SetCellsAmplitudeFraction(oc->GetCellsAmplitudeFraction());
284  dc->SetID(oc->GetID());
285  dc->SetDispersion(oc->GetDispersion());
286  dc->SetEmcCpvDistance(-1);
287  dc->SetChi2(-1);
288  dc->SetTOF(oc->GetTOF()); //time-of-flight
289  dc->SetNExMax(oc->GetNExMax()); //number of local maxima
290  dc->SetM02(oc->GetM02());
291  dc->SetM20(oc->GetM20());
292  dc->SetDistanceToBadChannel(oc->GetDistanceToBadChannel());
293  dc->SetMCEnergyFraction(oc->GetMCEnergyFraction());
294 
295  //MC
296  UInt_t nlabels = oc->GetNLabels();
297  Int_t *labels = oc->GetLabels();
298 
299  if (nlabels == 0 || !labels)
300  continue;
301 
302  AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
303  if (esdClus) {
304  TArrayI parents(nlabels, labels);
305  esdClus->AddLabels(parents);
306  }
307  else {
308  AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
309  if (aodClus)
310  aodClus->SetLabel(labels, nlabels);
311  }
312  }
313 }
314 
315 //________________________________________________________________________
317 {
318  // Clusterize
319 
320  if (fSubBackground) {
321  fClusterizer->SetInputCalibrated(kTRUE);
322  fClusterizer->SetCalibrationParameters(0);
323  }
324 
325  fClusterizer->Digits2Clusters("");
326 
327  if (fSubBackground) {
328  if (fCalibData) {
329  fClusterizer->SetInputCalibrated(kFALSE);
330  fClusterizer->SetCalibrationParameters(fCalibData);
331  }
332  }
333 }
334 
335 //________________________________________________________________________
337 {
338  // Fill digits array
339 
340  fDigitsArr->Clear("C");
341  switch (fInputCellType) {
342 
343  case kFEEData :
344  case kFEEDataMCOnly :
345  case kFEEDataExcludeMC :
346  {
347  // In case of MC productions done before aliroot tag v5-02-Rev09
348  // assing the cluster label to all the cells belonging to this cluster
349  // very rough
350  // Copied and simplified from AliEMCALTenderSupply
352  {
353  for (Int_t i = 0; i < fgkTotalCellNumber; i++)
354  {
355  fCellLabels [i] =-1 ;
356  fOrgClusterCellId[i] =-1 ;
357  }
358 
359  Int_t nClusters = InputEvent()->GetNumberOfCaloClusters();
360  for (Int_t i = 0; i < nClusters; i++)
361  {
362  AliVCluster *clus = InputEvent()->GetCaloCluster(i);
363 
364  if (!clus) continue;
365 
366  if (!clus->IsEMCAL()) continue ;
367 
368  Int_t label = clus->GetLabel();
369  UShort_t * index = clus->GetCellsAbsId() ;
370 
371  for(Int_t icell=0; icell < clus->GetNCells(); icell++)
372  {
374  fCellLabels[index[icell]] = label;
375 
376  fOrgClusterCellId[index[icell]] = i ; // index of the original cluster
377  } // cell in cluster loop
378  } // cluster loop
379  }
380 
381  Double_t avgE = 0; // for background subtraction
382  const Int_t ncells = fCaloCells->GetNumberOfCells();
383  for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell)
384  {
385  Float_t cellAmplitude=0;
386  Double_t cellTime=0, amp = 0, cellEFrac = 0;
387  Short_t cellNumber=0;
388  Int_t cellMCLabel=-1;
389  if (fCaloCells->GetCell(icell, cellNumber, amp, cellTime, cellMCLabel, cellEFrac) != kTRUE)
390  break;
391 
392  cellAmplitude = amp; // compilation problem
393 
394  if (fSetCellMCLabelFromCluster) cellMCLabel = fCellLabels[cellNumber];
395 
396  if (cellMCLabel > 0 && cellEFrac < 1e-6)
397  cellEFrac = 1;
398 
399  if (cellAmplitude < 1e-6 || cellNumber < 0)
400  continue;
401 
403  if (cellMCLabel <= 0)
404  continue;
405  else {
406  cellAmplitude *= cellEFrac;
407  cellEFrac = 1;
408  }
409  }
410  else if (fInputCellType == kFEEDataExcludeMC) {
411  if (cellMCLabel > 0)
412  continue;
413  else {
414  cellAmplitude *= 1 - cellEFrac;
415  cellEFrac = 0;
416  }
417  }
418 
419  if(!AcceptCell(cellNumber)) continue;
420 
421  //
422  // New way to set the cell MC labels,
423  // valid only for MC productions with aliroot > v5-07-21
424  //
425  TArrayI labeArr(0);
426  TArrayF eDepArr(0);
427  Int_t nLabels = 0;
428 
429  if(fSetCellMCLabelFromEdepFrac && fOrgClusterCellId[cellNumber] >= 0) // index can be negative if noisy cell that did not form cluster
430  {
431  cellMCLabel = -1;
432 
433  fCellLabels[cellNumber] = idigit;
434 
435  Int_t iclus = fOrgClusterCellId[cellNumber];
436 
437  if(iclus < 0)
438  {
439  AliInfo("Negative original cluster index, skip \n");
440  continue;
441  }
442 
443  AliVCluster *clus = InputEvent()->GetCaloCluster(iclus);
444 
445  for(Int_t icluscell=0; icluscell < clus->GetNCells(); icluscell++ )
446  {
447  if(cellNumber != clus->GetCellAbsId(icluscell)) continue ;
448 
449  // Select the MC label contributing, only if enough energy deposition
450 
451  fRecoUtils->RecalculateCellLabelsRemoveAddedGenerator(cellNumber, clus, MCEvent(),cellAmplitude, labeArr, eDepArr);
452  nLabels = labeArr.GetSize();
453  }
454  }
455 
456  AliEMCALDigit *digit = new((*fDigitsArr)[idigit]) AliEMCALDigit(cellMCLabel, cellMCLabel, cellNumber,
457  cellAmplitude, (Float_t)cellTime,
458  AliEMCALDigit::kHG,idigit, 0, 0, cellEFrac*cellAmplitude);
459 
460  if(nLabels > 0)
461  {
462  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
463  }
464 
465 
467  {
468  Float_t energy = cellAmplitude;
469  Float_t time = cellTime;
470  fClusterizer->Calibrate(energy,time,cellNumber);
471  digit->SetAmplitude(energy);
472  avgE += energy;
473  }
474 
475  idigit++;
476  }
477 
478  if (fSubBackground) {
479  avgE /= fGeom->GetNumberOfSuperModules()*48*24;
480  Int_t ndigis = fDigitsArr->GetEntries();
481  for (Int_t i = 0; i < ndigis; ++i) {
482  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
483  Double_t energy = digit->GetAmplitude() - avgE;
484  if (energy<=0.001) {
485  digit->SetAmplitude(0);
486  } else {
487  digit->SetAmplitude(energy);
488  }
489  }
490  }
491  }
492  break;
493 
494  case kPattern :
495  {
496  // Fill digits from a pattern
497  Int_t maxd = fGeom->GetNCells() / 4;
498  for (Int_t idigit = 0; idigit < maxd; idigit++){
499  if (idigit % 24 == 12) idigit += 12;
500  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
501  digit->SetId(idigit * 4);
502  digit->SetTime(600);
503  digit->SetTimeR(600);
504  digit->SetIndexInList(idigit);
505  digit->SetType(AliEMCALDigit::kHG);
506  digit->SetAmplitude(0.1);
507  }
508  }
509  break;
510 
511  case kL0FastORs :
512  case kL0FastORsTC :
513  case kL1FastORs :
514  {
515  // Fill digits from FastORs
516 
517  AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
518 
519  if (!triggers || !(triggers->GetEntries() > 0))
520  return;
521 
522  Int_t idigit = 0;
523  triggers->Reset();
524 
525  while ((triggers->Next())) {
526  Float_t L0Amplitude = 0;
527  triggers->GetAmplitude(L0Amplitude);
528 
529  if (L0Amplitude <= 0 && fInputCellType != kL1FastORs)
530  continue;
531 
532  Int_t L1Amplitude = 0;
533  triggers->GetL1TimeSum(L1Amplitude);
534 
535  if (L1Amplitude <= 0 && fInputCellType == kL1FastORs)
536  continue;
537 
538  Int_t triggerTime = 0;
539  Int_t ntimes = 0;
540  triggers->GetNL0Times(ntimes);
541 
542  if (ntimes < 1 && fInputCellType == kL0FastORsTC)
543  continue;
544 
545  if (ntimes > 0) {
546  Int_t trgtimes[25];
547  triggers->GetL0Times(trgtimes);
548  triggerTime = trgtimes[0];
549  }
550 
551  Int_t triggerCol = 0, triggerRow = 0;
552  triggers->GetPosition(triggerCol, triggerRow);
553 
554  Int_t find = -1;
555  fGeom->GetAbsFastORIndexFromPositionInEMCAL(triggerCol, triggerRow, find);
556 
557  if (find < 0)
558  continue;
559 
560  Int_t cidx[4] = {-1};
561  Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
562 
563  if (!ret)
564  continue;
565 
566  Float_t triggerAmplitude = 0;
567 
568  if (fInputCellType == kL1FastORs) {
569  triggerAmplitude = 0.25 * L1Amplitude; // it will add 4 cells for 1 amplitude
570  }
571  else {
572  triggerAmplitude = L0Amplitude; // 10 bit truncated, so it is already divided by 4
573  }
574 
575  for (Int_t idxpos = 0; idxpos < 4; idxpos++) {
576  Int_t triggerNumber = cidx[idxpos];
577  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
578  digit->SetId(triggerNumber);
579  digit->SetTime(triggerTime);
580  digit->SetTimeR(triggerTime);
581  digit->SetIndexInList(idigit);
582  digit->SetType(AliEMCALDigit::kHG);
583  digit->SetAmplitude(triggerAmplitude);
584  idigit++;
585  }
586  }
587  }
588  break;
589  }
590 }
591 
592 //________________________________________________________________________________________
594 
595  Bool_t accept = kTRUE;
596  if(fRejectExoticCells) {
597  //Remove exotic cells before making digits
598  Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
599  fRecoUtils->SwitchOnRejectExoticCell();//switch on and off
600  Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
601  Bool_t isEx = fRecoUtils->IsExoticCell(cellNumber, fCaloCells, bunchCrossNo);
602  if(isEx) accept = kFALSE;
603  if(!exRemoval) fRecoUtils->SwitchOffRejectExoticCell();//switch on and off
604  }
605  return accept;
606 }
607 
608 //________________________________________________________________________________________
610 {
611  // Go through clusters one by one and process separate correction
612  // as those were defined or not
613 
614  Int_t nclusters = fCaloClusters->GetEntriesFast();
615  for (Int_t icluster=0; icluster < nclusters; ++icluster) {
616  AliVCluster *clust = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
617  if (!clust)
618  continue;
619  if (!clust->IsEMCAL())
620  continue;
621 
622  // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
624  // careful, the the ClusterContainsBadChannel is dependent on
625  // SwitchOnBadChannelsRemoval, switching it ON automatically
626  // and returning to original value after processing
629 
630  Bool_t badResult = fRecoUtils->ClusterContainsBadChannel(fGeom, clust->GetCellsAbsId(), clust->GetNCells());
631 
632  // switch the bad channels removal back
633  if (!badRemoval)
635 
636  if (badResult) {
637  delete fCaloClusters->RemoveAt(icluster);
638  continue; //TODO is it really needed to remove it? Or should we flag it?
639  }
640  }
641 
642  // REMOVE EXOTIC CLUSTERS -------------------------------------
643  // does process local cell recalibration energy and time without replacing
644  // the global cell values, in case of no cell recalib done yet
645  if (fRejectExoticClusters) {
646  // careful, the IsExoticCluster is dependent on
647  // SwitchOnRejectExoticCell, switching it ON automatically
648  // and returning to original value after processing
649  Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
651 
652  // get bunch crossing
653  Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
654 
655  Bool_t exResult = fRecoUtils->IsExoticCluster(clust, fCaloCells, bunchCrossNo);
656 
657  // switch the exotic channels removal back
658  if (!exRemoval)
660 
661  if (exResult) {
662  delete fCaloClusters->RemoveAt(icluster);
663  continue; //TODO is it really needed to remove it? Or should we flag it?
664  }
665  }
666 
667  // FIDUCIAL CUT -----------------------------------------------
668  if (fFiducial) {
669  // depends on SetNumberOfCellsFromEMCALBorder
670  // SwitchOnNoFiducialBorderInEMCALEta0
672  delete fCaloClusters->RemoveAt(icluster);
673  continue; //TODO it would be nice to store the distance
674  }
675  }
676 
677  // NONLINEARITY -----------------------------------------------
678  if (fDoNonLinearity) {
679  Float_t correctedEnergy = fRecoUtils->CorrectClusterEnergyLinearity(clust);
680  clust->SetE(correctedEnergy);
681  }
682 
683  // DISTANCE TO BAD CHANNELS -----------------------------------
686  }
687 
688  fCaloClusters->Compress();
689 }
690 
691 //________________________________________________________________________________________
692 void AliAnalysisTaskEMCALClusterizeFast::TrackClusterMatching(AliVCluster *c, TClonesArray *tarr)
693 {
694  Float_t g[3]={0};
695  c->GetPosition(g);
696  TVector3 gpos(g);
697 
698  Double_t dEtaMin = 1e9;
699  Double_t dPhiMin = 1e9;
700  Int_t imin = -1;
701  Double_t ceta = gpos.Eta();
702  Double_t cphi = gpos.Phi();
703  const Int_t ntrks = tarr->GetEntries();
704  for(Int_t t = 0; t<ntrks; ++t) {
705  AliVTrack *track = static_cast<AliVTrack*>(tarr->At(t));
706  if (!track)
707  continue;
708  const AliExternalTrackParam *outp = track->GetOuterParam();
709  if (!outp)
710  continue;
711  Double_t trkPos[3] = {0.,0.,0.};
712  if (!outp->GetXYZ(trkPos))
713  continue;
714  TVector3 vec(trkPos[0],trkPos[1],trkPos[2]);
715  Double_t veta = vec.Eta();
716  Double_t vphi = vec.Phi();
717  if(vphi<0)
718  vphi += 2*TMath::Pi();
719  if (TMath::Abs(veta)>0.75 || (vphi<70*TMath::DegToRad()) || (vphi>190*TMath::DegToRad()))
720  continue;
721  Double_t dR = vec.DeltaR(gpos);
722  if(dR > 25)
723  continue;
724  Float_t tmpEta=0, tmpPhi=0;
725  if (0) {
726  AliExternalTrackParam trkParTemp(*outp); // retrieve the starting point every time before the extrapolation
727  Bool_t ret = fRecoUtils->ExtrapolateTrackToCluster(&trkParTemp, c, fRecoUtils->GetMass(), fRecoUtils->GetStep(), tmpEta, tmpPhi);
728  if (!ret)
729  continue;
730  } else {
731  tmpEta = ceta - veta;
732  tmpPhi = cphi - vphi;
733  }
734  if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMin) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMin)) {
735  dEtaMin = tmpEta;
736  dPhiMin = tmpPhi;
737  imin = t;
738  }
739  }
740  c->SetEmcCpvDistance(imin);
741  c->SetTrackDistance(dPhiMin, dEtaMin);
742 }
743 
744 //________________________________________________________________________________________
746 {
747  // Cluster energy, global position, cells and their amplitude fractions are restored.
748 
749  // tracks array for track/cluster matching
750  TClonesArray *tarr = 0;
751  if (!fTrackName.IsNull()) {
752  tarr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackName));
753  if (!tarr) {
754  AliError(Form("Cannot get tracks named %s", fTrackName.Data()));
755  }
756  }
757 
758  const Int_t Ncls = fClusterArr->GetEntries();
759  AliDebug(1, Form("total no of clusters %d", Ncls));
760 
761  for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i)
762  {
763  AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
764 
765  Int_t ncellsTrue = 0;
766  const Int_t ncells = recpoint->GetMultiplicity();
767  UShort_t absIds[ncells];
768  Double32_t ratios[ncells];
769  Int_t *dlist = recpoint->GetDigitsList();
770  Float_t *elist = recpoint->GetEnergiesList();
771  Double_t mcEnergy = 0;
772 
773  for (Int_t c = 0; c < ncells; ++c)
774  {
775  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
776  absIds[ncellsTrue] = digit->GetId();
777  ratios[ncellsTrue] = elist[c]/digit->GetAmplitude();
778 
779  if (digit->GetIparent(1) > 0)
780  mcEnergy += digit->GetDEParent(1)/recpoint->GetEnergy();
781 
782  ++ncellsTrue;
783  }
784 
785  if (ncellsTrue < 1)
786  {
787  AliWarning("Skipping cluster with no cells");
788  continue;
789  }
790 
791  // calculate new cluster position
792  TVector3 gpos;
793  recpoint->GetGlobalPosition(gpos);
794  Float_t g[3];
795  gpos.GetXYZ(g);
796 
797  AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
798 
799  AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
800  c->SetType(AliVCluster::kEMCALClusterv1);
801  c->SetE(recpoint->GetEnergy());
802  c->SetPosition(g);
803  c->SetNCells(ncellsTrue);
804  c->SetCellsAbsId(absIds);
805  c->SetCellsAmplitudeFraction(ratios);
806  c->SetID(recpoint->GetUniqueID());
807  c->SetDispersion(recpoint->GetDispersion());
808  c->SetEmcCpvDistance(-1);
809  c->SetChi2(-1);
810  c->SetTOF(recpoint->GetTime()) ; //time-of-flight
811  c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
812  Float_t elipAxis[2];
813  recpoint->GetElipsAxis(elipAxis);
814  c->SetM02(elipAxis[0]*elipAxis[0]);
815  c->SetM20(elipAxis[1]*elipAxis[1]);
816  c->SetMCEnergyFraction(mcEnergy);
817 
818  //
819  // MC labels
820  //
821  Int_t parentMult = 0;
822  Int_t *parentList = recpoint->GetParents(parentMult);
823  Float_t *parentListDE = recpoint->GetParentsDE(); // deposited energy
824 
825  c->SetLabel(parentList, parentMult);
826 
828  c->SetClusterMCEdepFractionFromEdepArray(parentListDE);
829 
830  //
831  // Set the cell energy deposition fraction map:
832  //
833  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
834  {
835  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
836 
837  // Get the digit that originated this cell cluster
838  //AliVCaloCells* cells = InputEvent()->GetEMCALCells();
839 
840  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
841  {
842  Int_t idigit = fCellLabels[absIds[icell]];
843 
844  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
845 
846  // Find the 4 MC labels that contributed to the cluster and their
847  // deposited energy in the current digit
848 
849  mcEdepFracPerCell[icell] = 0; // init
850 
851  Int_t nparents = dig->GetNiparent();
852  if ( nparents > 0 )
853  {
854  Int_t digLabel =-1 ;
855  Float_t edep = 0 ;
856  Float_t edepTot = 0 ;
857  Float_t mcEDepFrac[4] = {0,0,0,0};
858 
859  // all parents in digit
860  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
861  {
862  digLabel = dig->GetIparent (jndex+1);
863  edep = dig->GetDEParent(jndex+1);
864  edepTot += edep;
865 
866  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
867  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
868  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
869  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
870  } // all prarents in digit
871 
872  // Divide energy deposit by total deposited energy
873  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
874  if(edepTot > 0.01)
875  {
876  mcEdepFracPerCell[icell] = c->PackMCEdepFraction(mcEDepFrac);
877  }
878  } // at least one parent label in digit
879  } // cell in cluster loop
880 
881  c->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
882 
883  delete [] mcEdepFracPerCell;
884 
885  } // at least one parent in cluster, do the cell primary packing
886 
887  //
888  // Track matching
889  //
890  if (tarr)
891  TrackClusterMatching(c, tarr);
892  }
893 }
894 
895 //________________________________________________________________________
897 {
898  // Update cells in case re-calibration was done.
899  if (!fCalibData&&!fSubBackground)
900  return;
901 
902  const Int_t ncells = fCaloCells->GetNumberOfCells();
903  const Int_t ndigis = fDigitsArr->GetEntries();
904  if (ncells!=ndigis) {
905  fCaloCells->DeleteContainer();
906  fCaloCells->CreateContainer(ndigis);
907  }
908  for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
909  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
910  Double_t cellAmplitude = digit->GetCalibAmp();
911  Short_t cellNumber = digit->GetId();
912  Double_t cellTime = digit->GetTime();
913  fCaloCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
914  }
915 }
916 
917 //________________________________________________________________________
919 {
920  // Update cells in case re-calibration was done.
921 
922  const Int_t nents = fCaloClusters->GetEntries();
923  for (Int_t i=0;i<nents;++i) {
924  AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
925  if (!c)
926  continue;
927  if (c->IsEMCAL())
928  delete fCaloClusters->RemoveAt(i);
929  }
930 
931  fCaloClusters->Compress();
932 
934 }
935 
936 //________________________________________________________________________________________
938 {
939  // Select clusterization/unfolding algorithm and set all the needed parameters.
940 
941  if (InputEvent()->GetRunNumber()==fRun)
942  return;
943  fRun = InputEvent()->GetRunNumber();
944 
945  if (fJustUnfold){
946  // init the unfolding afterburner
947  delete fUnfolder;
948  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
949  return;
950  }
951 
952  if (fGeomName.Length()>0)
953  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
954  else
955  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
956  if (!fGeom) {
957  AliFatal("Geometry not available!!!");
958  return;
959  }
960 
961  if (!fGeomMatrixSet) {
962  if (fLoadGeomMatrices) {
963  for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
964  if (fGeomMatrix[mod]){
965  if (DebugLevel() > 2)
966  fGeomMatrix[mod]->Print();
967  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
968  }
969  }
970  } else { // get matrix from file (work around bug in aliroot)
971  for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
972  const TGeoHMatrix *gm = 0;
973  if (fEsd) {
974  gm = fEsd->GetEMCALMatrix(mod);
975  } else {
976  AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(fAod->GetHeader());
977  if(!aodheader) AliFatal("Not a standard AOD");
978  if (aodheader) {
979  gm = aodheader->GetEMCALMatrix(mod);
980  }
981  }
982  if (gm) {
983  if (DebugLevel() > 2)
984  gm->Print();
985  fGeom->SetMisalMatrix(gm,mod);
986  }
987  }
988  }
989  fGeomMatrixSet=kTRUE;
990  }
991 
992  // setup digit array if needed
993  if (!fDigitsArr) {
994  fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
995  fDigitsArr->SetOwner(1);
996  }
997 
998  // then setup clusterizer
999  if (fClusterizer) {
1000  // avoid to delete digits array
1001  fClusterizer->SetDigitsArr(0);
1002  delete fClusterizer;
1003  }
1004  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1005  fClusterizer = new AliEMCALClusterizerv1(fGeom);
1006  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
1007  AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
1008  clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
1009  clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
1010  fClusterizer = clusterizer;
1011  }
1012  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1013  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1014  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv3)
1015  fClusterizer = new AliEMCALClusterizerv3(fGeom);
1016  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
1017  AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
1018  clusterizer->SetNphi(fNPhi);
1019  clusterizer->SetNeta(fNEta);
1020  clusterizer->SetShiftPhi(fShiftPhi);
1021  clusterizer->SetShiftEta(fShiftEta);
1022  clusterizer->SetTRUshift(fTRUShift);
1023  fClusterizer = clusterizer;
1024  }
1025  else {
1026  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1027  }
1028  fClusterizer->InitParameters(fRecParam);
1029 
1030  if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
1031  AliCDBManager *cdb = AliCDBManager::Instance();
1032  if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
1033  cdb->SetDefaultStorage(fOCDBpath);
1034  if (fRun!=cdb->GetRun())
1035  cdb->SetRun(fRun);
1036  }
1037  if (!fCalibData&&fLoadCalib&&fRun>0) {
1038  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
1039  if (entry)
1040  fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
1041  if (!fCalibData)
1042  AliFatal("Calibration parameters not found in CDB!");
1043  }
1044  if (!fPedestalData&&fLoadPed&&fRun>0) {
1045  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
1046  if (entry)
1047  fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
1048  }
1049  if (fCalibData) {
1050  fClusterizer->SetInputCalibrated(kFALSE);
1051  fClusterizer->SetCalibrationParameters(fCalibData);
1052  } else {
1053  fClusterizer->SetInputCalibrated(kTRUE);
1054  }
1055  fClusterizer->SetCaloCalibPedestal(fPedestalData);
1056  fClusterizer->SetJustClusters(kTRUE);
1057  fClusterizer->SetDigitsArr(fDigitsArr);
1058  fClusterizer->SetOutput(0);
1059  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1060 
1061  // Get the emcal cells
1063  if (fCaloCellsName.IsNull()) {
1064  fCaloCells = InputEvent()->GetEMCALCells();
1065  }
1066  else {
1067  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
1068  if (!fCaloCells)
1069  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
1070  }
1071  if (!fCaloCells)
1072  AliFatal("Could not get EMCal cells!");
1073  }
1074 
1075  // Set output clusters collection
1076  if (!fAttachClusters) {
1078  return;
1079  }
1080 
1081  if (!fCaloClusters) {
1082  if (fCaloClustersName.IsNull()) { //overwrite mode
1083  if (fEsd)
1084  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
1085  else if (fAod)
1086  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
1087  }
1088  else {
1089  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloClustersName));
1090 
1091  if (!fCaloClusters) {
1092  if (fEsd)
1093  fCaloClusters = new TClonesArray("AliESDCaloCluster");
1094  else if (fAod)
1095  fCaloClusters = new TClonesArray("AliAODCaloCluster");
1096 
1097  fCaloClusters->SetName(fCaloClustersName);
1098  InputEvent()->AddObject(fCaloClusters);
1099  }
1100  }
1101 
1102  if (!fCaloClusters)
1103  AliFatal("Could not get cluster collection!");
1104 
1105  TClass *cl = fCaloClusters->GetClass();
1106  if (!cl->GetBaseClass("AliVCluster")) {
1107  AliFatal(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloClusters->GetName()));
1108  fCaloClusters = 0;
1109  return;
1110  }
1111  }
1112 }
void RecalculateCellLabelsRemoveAddedGenerator(Int_t absID, AliVCluster *clus, AliMCEvent *mc, Float_t &amp, TArrayI &labeArr, TArrayF &eDepArr) const
void Print(std::ostream &o, const char *name, Double_t dT, Double_t dVM, Double_t alldT, Double_t alldVM)
Definition: PlotSysInfo.C:121
double Double_t
Definition: External.C:58
TString fOutputAODBrName
AOD Branch with output clusters.
Bool_t CheckCellFiducialRegion(const AliEMCALGeometry *geom, const AliVCluster *cluster, AliVCaloCells *cells)
Double_t GetMass() const
energy
Definition: HFPtSpectrum.C:44
Bool_t IsExoticCluster(const AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=0)
TCanvas * c
Definition: TestFitELoss.C:172
AliEMCALAfterBurnerUF * fUnfolder
clusterizer
Bool_t IsBadChannelsRemovalSwitchedOn() const
Bool_t IsExoticCell(Int_t absId, AliVCaloCells *cells, Int_t bc=-1)
void RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
virtual void TrackClusterMatching(AliVCluster *c, TClonesArray *tarr)
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
void SwitchOffBadChannelsRemoval()
virtual void CopyClusters(TClonesArray *orig, TClonesArray *dest)
short Short_t
Definition: External.C:23
static Bool_t ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, const AliVCluster *cluster, Double_t mass, Double_t step, Float_t &tmpEta, Float_t &tmpPhi)
Double_t GetStep() const
TGeoHMatrix * fGeomMatrix[AliEMCALGeoParams::fgkEMCALModules]
Bool_t IsRejectExoticCell() const
unsigned short UShort_t
Definition: External.C:28
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Float_t CorrectClusterEnergyLinearity(AliVCluster *clu)
Bool_t ClusterContainsBadChannel(const AliEMCALGeometry *geom, const UShort_t *cellList, Int_t nCells)