AliPhysics  master (3d17d9d)
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 
827  c->SetClusterMCEdepFractionFromEdepArray(parentListDE);
828 
829  //
830  // Set the cell energy deposition fraction map:
831  //
832  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
833  {
834  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
835 
836  // Get the digit that originated this cell cluster
837  //AliVCaloCells* cells = InputEvent()->GetEMCALCells();
838 
839  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
840  {
841  Int_t idigit = fCellLabels[absIds[icell]];
842 
843  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
844 
845  // Find the 4 MC labels that contributed to the cluster and their
846  // deposited energy in the current digit
847 
848  mcEdepFracPerCell[icell] = 0; // init
849 
850  Int_t nparents = dig->GetNiparent();
851  if ( nparents > 0 )
852  {
853  Int_t digLabel =-1 ;
854  Float_t edep = 0 ;
855  Float_t edepTot = 0 ;
856  Float_t mcEDepFrac[4] = {0,0,0,0};
857 
858  // all parents in digit
859  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
860  {
861  digLabel = dig->GetIparent (jndex+1);
862  edep = dig->GetDEParent(jndex+1);
863  edepTot += edep;
864 
865  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
866  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
867  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
868  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
869  } // all prarents in digit
870 
871  // Divide energy deposit by total deposited energy
872  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
873  if(edepTot > 0.01)
874  {
875  mcEdepFracPerCell[icell] = c->PackMCEdepFraction(mcEDepFrac);
876  }
877  } // at least one parent label in digit
878  } // cell in cluster loop
879 
880  c->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
881 
882  delete [] mcEdepFracPerCell;
883 
884  } // at least one parent in cluster, do the cell primary packing
885 
886  //
887  // Track matching
888  //
889  if (tarr)
890  TrackClusterMatching(c, tarr);
891  }
892 }
893 
894 //________________________________________________________________________
896 {
897  // Update cells in case re-calibration was done.
898  if (!fCalibData&&!fSubBackground)
899  return;
900 
901  const Int_t ncells = fCaloCells->GetNumberOfCells();
902  const Int_t ndigis = fDigitsArr->GetEntries();
903  if (ncells!=ndigis) {
904  fCaloCells->DeleteContainer();
905  fCaloCells->CreateContainer(ndigis);
906  }
907  for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
908  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
909  Double_t cellAmplitude = digit->GetCalibAmp();
910  Short_t cellNumber = digit->GetId();
911  Double_t cellTime = digit->GetTime();
912  fCaloCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
913  }
914 }
915 
916 //________________________________________________________________________
918 {
919  // Update cells in case re-calibration was done.
920 
921  const Int_t nents = fCaloClusters->GetEntries();
922  for (Int_t i=0;i<nents;++i) {
923  AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
924  if (!c)
925  continue;
926  if (c->IsEMCAL())
927  delete fCaloClusters->RemoveAt(i);
928  }
929 
930  fCaloClusters->Compress();
931 
933 }
934 
935 //________________________________________________________________________________________
937 {
938  // Select clusterization/unfolding algorithm and set all the needed parameters.
939 
940  if (InputEvent()->GetRunNumber()==fRun)
941  return;
942  fRun = InputEvent()->GetRunNumber();
943 
944  if (fJustUnfold){
945  // init the unfolding afterburner
946  delete fUnfolder;
947  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
948  return;
949  }
950 
951  if (fGeomName.Length()>0)
952  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
953  else
954  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
955  if (!fGeom) {
956  AliFatal("Geometry not available!!!");
957  return;
958  }
959 
960  if (!fGeomMatrixSet) {
961  if (fLoadGeomMatrices) {
962  for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
963  if (fGeomMatrix[mod]){
964  if (DebugLevel() > 2)
965  fGeomMatrix[mod]->Print();
966  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
967  }
968  }
969  } else { // get matrix from file (work around bug in aliroot)
970  for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
971  const TGeoHMatrix *gm = 0;
972  if (fEsd) {
973  gm = fEsd->GetEMCALMatrix(mod);
974  } else {
975  AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(fAod->GetHeader());
976  if(!aodheader) AliFatal("Not a standard AOD");
977  if (aodheader) {
978  gm = aodheader->GetEMCALMatrix(mod);
979  }
980  }
981  if (gm) {
982  if (DebugLevel() > 2)
983  gm->Print();
984  fGeom->SetMisalMatrix(gm,mod);
985  }
986  }
987  }
988  fGeomMatrixSet=kTRUE;
989  }
990 
991  // setup digit array if needed
992  if (!fDigitsArr) {
993  fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
994  fDigitsArr->SetOwner(1);
995  }
996 
997  // then setup clusterizer
998  if (fClusterizer) {
999  // avoid to delete digits array
1000  fClusterizer->SetDigitsArr(0);
1001  delete fClusterizer;
1002  }
1003  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1004  fClusterizer = new AliEMCALClusterizerv1(fGeom);
1005  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
1006  AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
1007  clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
1008  clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
1009  fClusterizer = clusterizer;
1010  }
1011  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1012  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1013  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv3)
1014  fClusterizer = new AliEMCALClusterizerv3(fGeom);
1015  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
1016  AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
1017  clusterizer->SetNphi(fNPhi);
1018  clusterizer->SetNeta(fNEta);
1019  clusterizer->SetShiftPhi(fShiftPhi);
1020  clusterizer->SetShiftEta(fShiftEta);
1021  clusterizer->SetTRUshift(fTRUShift);
1022  fClusterizer = clusterizer;
1023  }
1024  else {
1025  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1026  }
1027  fClusterizer->InitParameters(fRecParam);
1028 
1029  if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
1030  AliCDBManager *cdb = AliCDBManager::Instance();
1031  if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
1032  cdb->SetDefaultStorage(fOCDBpath);
1033  if (fRun!=cdb->GetRun())
1034  cdb->SetRun(fRun);
1035  }
1036  if (!fCalibData&&fLoadCalib&&fRun>0) {
1037  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
1038  if (entry)
1039  fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
1040  if (!fCalibData)
1041  AliFatal("Calibration parameters not found in CDB!");
1042  }
1043  if (!fPedestalData&&fLoadPed&&fRun>0) {
1044  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
1045  if (entry)
1046  fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
1047  }
1048  if (fCalibData) {
1049  fClusterizer->SetInputCalibrated(kFALSE);
1050  fClusterizer->SetCalibrationParameters(fCalibData);
1051  } else {
1052  fClusterizer->SetInputCalibrated(kTRUE);
1053  }
1054  fClusterizer->SetCaloCalibPedestal(fPedestalData);
1055  fClusterizer->SetJustClusters(kTRUE);
1056  fClusterizer->SetDigitsArr(fDigitsArr);
1057  fClusterizer->SetOutput(0);
1058  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1059 
1060  // Get the emcal cells
1062  if (fCaloCellsName.IsNull()) {
1063  fCaloCells = InputEvent()->GetEMCALCells();
1064  }
1065  else {
1066  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
1067  if (!fCaloCells)
1068  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
1069  }
1070  if (!fCaloCells)
1071  AliFatal("Could not get EMCal cells!");
1072  }
1073 
1074  // Set output clusters collection
1075  if (!fAttachClusters) {
1077  return;
1078  }
1079 
1080  if (!fCaloClusters) {
1081  if (fCaloClustersName.IsNull()) { //overwrite mode
1082  if (fEsd)
1083  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
1084  else if (fAod)
1085  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
1086  }
1087  else {
1088  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloClustersName));
1089 
1090  if (!fCaloClusters) {
1091  if (fEsd)
1092  fCaloClusters = new TClonesArray("AliESDCaloCluster");
1093  else if (fAod)
1094  fCaloClusters = new TClonesArray("AliAODCaloCluster");
1095 
1096  fCaloClusters->SetName(fCaloClustersName);
1097  InputEvent()->AddObject(fCaloClusters);
1098  }
1099  }
1100 
1101  if (!fCaloClusters)
1102  AliFatal("Could not get cluster collection!");
1103 
1104  TClass *cl = fCaloClusters->GetClass();
1105  if (!cl->GetBaseClass("AliVCluster")) {
1106  AliFatal(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloClusters->GetName()));
1107  fCaloClusters = 0;
1108  return;
1109  }
1110  }
1111 }
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:45
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)